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A  detailed  understanding  of  ceramics  requires  constitutive  models 
for  grains  and  grain  boundaries  and  a  method  for  predicting  the  response 
of  ceramic  specimens  to  experimental  loads.  Voronoi  diagrams  are  used  to 
construct  realistic  grain  topologies,  and  continuum  damage  relations  are 
used  to  represent  the  behavior  of  grain  boundaries.  A  robust  solution 
algorithm  is  used  to  obtain  solutions  that  include  softening  or  failure. 
These  results  shov  the  evolution  of  microcracks,  crack  branching,  crack 
bridging,  and  the  ultimate  development  of  microcracks.  Although  the 
computational  domain  only  contains  a  small  number  of  grains,  the  results 
indicate  that  such  simulations  are  feasible  and  that  extension  of  the  anal¬ 
ysis  to  three  dimensions  is  a  logical  step. 
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1.  INTRODUCTION 

For  quasi-brittle  materials  a  significant  amount  of  research  has  been  conducted 
wider  the  assumption  that  a  process  zone  develops  ahead  of  a  crack  tip.  Here,  a  crack 
process  zone  is  defined  to  be  a  zone  of  distributed  microcracks  developing  prior  to  a 
macrocrack  in  a  region  both  above  and  below  the  plane  of  the  macrocrack.  The  very 
existence  of  such  a  zone  has  been  the  subject  of  major  inquiry  [Mindess,  1990]. 
Experimental  evidence  [Wang,  Schreyer  and  Rutland,  1990;  van  Mier,  1990]  indicates  that 
such  a  zone  is  extremely  hard  to  identify  or  is  nonexistent  A  theoretical  investigation 
[Horii  and  Niimalendran,  1990]  shows  that  even  if  a  zone  of  microcracks  develops  prior  to 
the  formation  of  a  macrocrack,  the  influence  of  the  zone  is  a  relatively  insignificant  part  of 
the  response  of  the  structure.  Load-deflection  curves  obtained  in  the  laboratory  and  post- 
limit  sectioning  of  samples  indicate  that  there  is  a  significant  amount  of  activity  prior  to  die 
formation  of  a  macrocrack  but  not  in  the  form  of  a  crack  process  zone.  Instead, 
microcracks  open  prior  to  the  macrocrack  with  the  result  that  certain  elements  such  as 
inclusions  in  ceramics  maintain  contact  longer  than  the  surrounding  matrix.  This  is  called 
the  crack  bridging  effect  or  "wake"  effect  which  is  interpreted  as  shielding  the  crack  tip 
from  applied  stress  [Reichl  and  Steinbrech,  1988].  Furthermore,  aggregates  or  inclusions 
cause  the  resulting  crack  surface  to  be  strongly  nonplanar  or  "tortuous,"  also  called  crack 
deflection  [Faber  and  Evans,  1983].  The  development  of  crack  bridging  [Kobayashi  and 
Shochey,  1987;  Li,  1990],  tortuosity  and  a  nonlinear  crack  front  formed  by  the  presence  of 
inclusions  causes  a  specimen  to  appear  less  brittle  than  an  identical  geometrical  specimen 
composed  of  tire  parent  matrix  material.  Such  materials  are  called  quasi -brittle  because 
classical  linear  elastic  fracture  mechanics  cannot  be  applied  without  modifications  to 
account  for  the  apparent  ductility  which  is  developed  prior  to  failure. 

The  most  straightforward  experimental  examples  involve  Mode  I  cracks  as 
developed,  for  example,  in  tension  or  three-point  and  four-point  bending  tests.  In  addition 
to  the  apparent  ductility  displayed  prior  to  the  peak  load,  a  softening  branch  is  exhibited  in 
the  post-peak  regime  of  the  load-deflection  curve  if  the  specimen  is  small  enough.  This 
softening  is  another  manifestation  of  apparent  ductility.  Under  certain  combinations  of 
material  and  geometrical  properties,  cracks  will  propagate  in  a  continuous  but  slow  rate 
under  displacement  control  [Maji  et  al.,  1990].  This  suggests  that  the  crack  driving  "force" 
is  just  sufficient  to  overcome  the  crack  "resistance"  and  that  the  process  is  not  necessarily 
unstable  as  is  commonly  claimed.  Suppose  that  cracks  may  grow  at  a  rate  governed  by  an 
inherent  material  property  and  the  crack  driving  force.  Then  if  a  specimen  is  loaded 
dynamically,  the  average  stress  could  significantly  exceed  tire  static  stress  necessary  to 
cause  failure,  and  any  number  of  pre-existing  flaws  could  serve  to  initiate  cracks. 
However,  most  dynamic  loads  are  impulsive  in  nature  so  the  growth  of  these  cracks  is 
soon  arrested  because  of  the  lack  of  a  driving  farce.  The  result  is  the  potential  existence  of 
a  region  of  microcracks.  Therefore,  for  dynamic  situations,  tire  evolution  of  a  macrocrack 
must  be  understood  when  considering  a  microcrack  process  zone.  For  quasi-static 
processes,  which  are  the  subject  of  this  research,  there  is  little  experimental  evidence  that 
displays  a  crack  process  zone  in  the  form  of  a  region  of  distributed  microcracks  prior  to  the 
macrocrack.  The  subject  is  one  of  considerable  controversy. 


2.  RESEARCH  OBJECTIVE 

Existing  experimental  results  indicate  that  crack  branching  occurs,  and  existing 
cracks  sustain  a  significant  amount  of  load  carrying  capacity.  This  feature  has  been 
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demonstrated  in  an  alumina  ceramic:  "Active  grain  bridges  were  observed  along  the  entire 
crack  trace  and  over  the  entire  propagation  distance.  No  indication  of  a  microcrack -cloud 
zone  was  observed.  ..."[Rodel.  1990].  To  capture  the  feature  of  support  along  the  entire 
crack  trace,  the  damage  crack  concept  is  proposed  in  this  research  in  which  the  traction  on 
the  crack  face  is  related  to  the  discontinuity  in  displacement  through  a  crack  modulus. 
Damage,  and  hence  the  modulus,  also  vary  with  the  discontinuity  in  displacement  across 
the  crack.  The  model  used  in  this  project  yielded  a  problem  somewhat  analogous  to  a 
spring  foundation  but  with  a  nonlinear  spring.  From  conventional  damage  concepts,  the 
stiffness  coefficients  were  related  to  the  amount  of  ligament  area  still  intact  ui  comparison  to 
the  total  area  of  die  crack.  As  the  crack  evolved,  the  stiffness  of  the  complete  specimen  was 
correlated  with  the  degree  of  damage  in  the  crack.  There  was  no  assumption  that 
distributed  microcracks  exist  with  this  approach.  The  predictions  of  stiffness  for  die 
specimen,  or  a  part  of  the  specimen,  and  the  strength  of  the  specimen  were  correlated  with 
experimental  data  in  the  literature  to  provide  an  initial  indication  of  the  magnitude  of  the 
crack  stiffness. 

Numerical  solutions  were  obtained  for  the  case  of  a  compact  tension  specimen.  The 
interface  was  modeled  as  a  damage  crack  with  a  small  amount  of  initial  damage.  The 
change  in  response  of  the  specimen  due  to  the  crack  propagation  provided  an  indirect 
measure  of  the  original  strength  of  the  bond.  Problems  involving  softening  and  localization 
are  notoriously  difficult  to  solve  because  the  tangent  stiffness  matrix  becomes  singular  at 
the  peak  load,  and  multiple  solution  paths  exist  for  displacement  beyond  those  at  peak  load. 
Existing  numerical  algorithms  were  attempted,  but  even  the  simplest  situation  of  a  two- 
dimensional  mode  for  uniaxial  tension  required  an  inordinate  amount  of  computer  time. 
One  of  the  more  important  achievements  of  this  research  was  the  development  of  a  robust 
and  efficient  algorithm  for  solving  problems  involving  softening  and  localization. 

The  issue  of  crack  branching  was  addressed  by  introducing  inclusions  to  cause 
multiple  interruptions.  The  reason  why  a  microcrack  ceases  to  grow  in  preference  to  an 
alternative  path  was  addressed  numerically  by  investigating  alternative  arrangements  of 
inclusions.  Specifically,  angular  shape  inclusions  were  considered  and  an  attempt  was 
made  to  show  how  inclusions  can  rotate  to  cause  bridging  by  interrupting  the  crack  paths. 

The  following  are  objectives  of  the  numerical  phase  of  this  study: 

1.  By  an  indirect  process  of  comparing  experimental  data  on  tensile  failure  with  predictions 
based  on  a  spectrum  of  grain  boundary  material  parameters,  we  will  be  able  to  determine 
die  mechanical  properties  of  the  grain  boundary  itself. 

2.  We  will  be  able  to  show  how  crack  branching  occurs.  These  results  will  be  important  for 
constructing  nonlocal  continuum  models  for  representing  die  ceramic. 

3.  We  will  be  able  to  demonstrate  that  a  micros tmctural  analysis  can  be  performed  so  that 
the  effects  of  additional  features  such  as  porosity  and  anisotropy  can  be  incorporated. 

4.  The  basis  will  be  established  for  looking  at  toughening  mechanisms  provided  by  fibers 
and  particles  in  a  systematic  manner  rather  than  by  empirical  methods  currently  in  use 
[Ptzzotti  et  al.  1990a,  1990b]. 

5.  The  ultimate  advantage  of  the  numerical  approach  will  be  the  capability  for  identifying 
quantitatively  those  individual  features  that  will  have  the  greatest  potential  for  enhancing 
strength  and  ductility. 

It  would  be  useful  for  analysts  if  a  ceramic  could  be  modeled  as  a  continuum, 
process  zone  mechanism  varying  from  one  ceramic  to  the  other,  as  exhibited,  for  example, 
by  toughening  in  which  fracture  resistance  systematically  increases  with  crack  extension. 
Such  features  appear  if  mechanisms  such  as  phase  transformations,  grain  bridging  and 
whisker  reinforcements  are  present  [Evans,  1990].  These  features  can  be  captured  in  a 
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continuum  model  only  if  softening  and  nonlocal  aspects  are  included.  The  nonlocal  pan  of 
a  constitutive  equation  controls  tne  size  of  the  predicted  process  zone  and  maintains  the 
mathematical  well-posedness  of  the  governing  equilibrium  equation,  i.e.,  the  original 
elliptic  nature  of  the  equation  remains  elliptic  as  a  crack  process  zone  is  established 

The  compact-tension  tests  of  alumina  ceramic  were  performed  on  large  specimens 
using  the  Instron  Testing  Machine,  and  on  a  much  smaller  scale  inside  the  chamber  of  a 
Scanning  Electron  Microscope  (SEM).  The  small  in-si tu  tests  involving  the  SEM  provided 
the  phenomenon  of  microfracturing  both  for  scientific  understanding  and  for  practical 
application.  In  this  case,  the  special  loading  stage,  which  can  be  mounted  in  die  SEM 
chamber,  was  designed  including  the  closed-loop  feedback  to  the  loading  device. 

In  summary,  die  objectives  of  die  experimental  phase  were: 

(i)  to  determine  the  mechanics  of  bridging 

(ii)  to  provide  a  real-time  description  of  macrocrack  interaction  with  microcracks 

(iii)  to  provide  those  measurements  useful  for  verifying  the  theoretical  approach  such  as 
cracking-opening  displacements,  loads,  displacements  at  the  point  of  application  of  the 
load,  crack  lengths  and  locations  of  cracks 

(iv)  to  continue  the  development  of  the  in-situ  testing  device  to  include  a  displacement 
closed-loop  control  feature. 


3.  SUMMARY  OF  RESEARCH 
3.1  Experimental 

It  is  generally  recognized  that  the  nature  of  grain  boundaries  and  microstructures 
affects  the  properties  and  behavior  of  ceramics.  It  is  important  to  understand  the  ceramic 
processing  characteristics  which  affect  the  nature  of  grain  boundaries.  Different  processing 
techniques  have  been  attempted  to  produce  the  desired  microstructure  of  alumina  in  the 
laboratory  (See  progress  report,  Nov.  1, 1992).  Because  the  size  of  a  ceramic  sample 
cannot  be  augmented  to  perform  mechanical  tests,  and  because  the  ceramic  is  not  fully 
dense,  an  alumina  bar  measuring  30  mm  X  30  mm  X  230  mm  was  obtained  from  the 
Coots  Technical  Ceramics  Co.  The  relative  density  of  the  alumina  is  99.5%.  (see 
Appendix  1  for  the  processing  and  properties  of  alumina). 

Ceramic  compact  tension  tests  were  conducted  on  an  Instron.  Since  the  stiffness 
difference  between  the  load  frame  of  Instron  and  the  compact  tension  specimen  is  large,  A 
reasonable  load-deflection  curve  could  not  be  obtained.  The  Scanning  Election 
Microscope  (SEM)  is  widely  used  as  an  analytical  tool  in  the  field  of  failure  and  fracture 
analysis,  primarily  because  of  its  combination  of  good  field  depth  and  high  resolution. 
For  the  studying  of  the  micofracture  process  of  the  ceramic,  the  closed-loop  control  in- 
situ  loading  stage,  which  can  be  mounted  inside  the  SEM  chamber,  was  designed  to 
perform  a  real  time  fracture  test 

Compared  with  most  metallic  materials,  the  mechanical  properties  of  engineering 
ceramics  are  much  more  sensitive  to  factors  such  as  size,  shape  and  surface  finish.  When 
tests  are  conducted,  machining  of  some  of  the  surface  is  required  to  remove  surface  flaws. 
Ceramic  materials  are  difficult  and  expensive  to  machine  due  to  their  high  hardness  and 
brittle  nature.  The  tool  must  have  a  higher  hardness  than  the  ceramic  being  machined,  and 
must  be  of  a  configuration  that  removes  surface  stock  without  overstressing  the 
component  For  these  reasons,  several  special  machines  were  purchased  to  prepare  the 
specimen. 
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Examination  of  fracture  surfaces  with  SEM  can  provide  useful  information  on  such 
aspects  as  local  fracture  mechanisms  and  fracture  propagation  direction.  After  the  test,  the 
fracture  surface  is  examined  by  the  SEM.  Since  ceramic  samples  lack  sufficient 
conductivity  for  analysis  in  SEM,  to  produce  a  conductive  layer,  carbon,  gold,  or 
aluminum  is  evaporated  onto  the  specimen  surface.  Fractographies  show  the  fracture  of 
alumina  is  intergranular.  Intergranular  fracture  is  usually  the  most  easily  recognizable 
mode  of  fracture,  but  determining  the  cause  may  be  quite  difficult 


3.2  Theoretical 

Inclusions  in  ceramics  provide  a  toughening  effect  by  rotating  and  impeding  crack 
openings  through  what  is  called  "crack  bridging."  A  major  objective  of  this  research  is  to 
investigate  in  detail  the  mechanisms  that  occur  in  track  bridging;  however,  a  study  of 
micrographs  available  in  the  literature  indicate  that  the  same  phenomenon  appears  with 
grains,  only  at  a  smaller  scale.  In  particular,  crack  branching  occurs  around  grains  with 
one  of  these  microcracks,  ultimately  evolving  to  a  marocrack,  with  large  grains  appearing 
as  "inclusions"  which  rotate  as  rigid  bodies.  Based  on  this  observation,  the  theoretical 
phase  involved  an  investigation  of  fundamental  properties  of  grain  boundaries.  Such  an 
analysis  can  be  closely  guided  by  in-situ  testing  in  a  scanning  electron  microscope  (SEM), 
and  theoretical  results  can  provide  insight  into  how  the  mechanical  properties  of  ceramics 
can  be  improved. 

Continuous  damage  mechanics  is  concerned  only  with  the  description  of 
progressive  weakening  of  solids  due  to  the  development  of  microcracks  and  microvoids; 
therefore,  continuum  damage  mechanics  is  die  most  appropriate  approach  for  analyzing  the 
nonlinear  behavior  of  quasi-brittle  materials.  Based  on  microstructural  analysis,  it  was 
proposed  to  numerically  model  in  two  dimensions  a  small  region  composed  of  a  set  of 
grains  and  grain  boundaries,  as  shown  in  Fig.  1.  Each  grain  was  modeled  as  an  elastic 
material  while  the  grain  boundary  was  considered  to  be  a  thin  region  of  elastic  damaging 
material.  An  initial  assumption  was  that  each  grain  is  scalar  isotropic  but  at  a  later  stage,  a 
two-parameter  isotropic  model  was  introduced. 

In  our  finite  element  codes,  dynamic  relaxation  (DR),  an  explicit  incremental 
iterative  method,  was  used  to  obtain  preliminary  numerical  results.  The  DR  method  is 
based  on  the  fact  that  the  static  solution  is  the  steady  state  part  erf  the  transient  response  for 
a  temporal-step  load  [Underwood,  1983].  This  method  is  especially  attractive  for  problems 
with  highly  nonlinear  geometric  and  material  behaviors,  which  include  limit  points  and 
regions  of  very  soft  stiffness  characteristics,  like  the  issue  in  our  research.  Crack 
nucleation  and  propagation  along  the  grain  boundary  are  shown  in  the  numerical  results. 
The  effects  of  microcracks,  voids  and  inclusions  are  also  observed. 

Numerical  investigations  are  difficult  because  of  the  ill -conditioning  inherent  in  the 
large  difference  between  the  characteristic  dimension  of  the  grain  and  the  thickness  of  the 
grain  boundary.  One  way  to  reduce  the  condition  number  is  to  model  each  grain  with  finite 
elements  of  dimension  similar  to  the  thickness  of  the  grain  boundary.  Tire  result  would  be 
an  inordinate  number  of  elements  that  could  not  be  handled  by  the  work  stations  available 
for  this  research.  Another  approach  is  to  invoke  constraints  to  remove  the  largest 
eigenvalues  which,  in  turn,  would  reduce  tire  condition  number  and  allow  for  a  possible 
numerical  solution  for  a  sufficiently  large  domain  so  that  meaningful  results  could  be 
obtained  [Hueck  and  Schreyer,  1992].  A  method  for  efficiently  incorporating  constraints 
resulted  in  a  paper  submitted  for  publication  [Schreyer  and  Parsons,  1994], 
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4.  DETAILS  OF  RESEARCH 


4.1  Ceramic  Properties  (Alumina) 

Ceramic  materials  include  oxides,  carbides,  sulfides,  and  intennetallic  compounds, 
which  are  joined  either  by  covalent  or  ionic  bonds.  Ceramic  moduli  are  generally  much 
larger  than  those  of  metals,  reflecting  greater  stiffness  of  bonds.  And  since  ceramics  are 
largely  composed  of  light  atoms  (oxygen,  carbon,  aluminum)  and  their  structures  are 
frequently  not  closely-packed,  their  densities  are  low.  Most  ceramics  are  crystalline,  but 
unlike  metals,  they  do  not  have  closely-packed  planes  on  which  dislocation  motion  can 
occur.  Therefore,  ceramic  materials  trad  to  be  very  brittle  compared  to  metals.  Typical 

ceramics  have  very  high  melting  temperatures  (like  alumina,  2050°C),  which  explains  their 
good  creep  properties.  Also,  many  these  materials  have  superior  wear  resistance,  and 
have  been  used  for  bearings  and  machine  tools.  Most  ceramics,  however,  are  too  brittle  for 
critical  loading-bearing  applications.  Consequently,  a  vast  amount  of  research  has  been 
devoted  to  improving  the  toughness  of  ceramics. 

Most  traditional  ceramics  are  monolithic  (single  phase)  and  have  very  low  fracture 
toughness.  Our  research  was  concentrated  on  monolithic  ceramics,  in  particular  aluminum 
oxide  (i.e.  alumina),  which  is  widely  used  in  industry. 

4.2  Microstructural  study  (alumina) 

Ceramics  are  formed  by  the  application  of  a  powder  to  high  temperatures  perhaps 
with  the  assistance  of  pressure  and  additives.  The  particles  sinter  together  into  a 
microstructure  consisting  of  individual  crystals  or  polycrystals,  i.e.,  grains  (made  up  of 
millions  of  small  crystals  [Ashby,  1986])  separated  by  grain  boundaries  and  residual 
porosity  (see  Fig.  2).  Each  grain  is  more  or  less  a  perfect  crystal.  The  grain  boundary  is 
obviously  complicated,  and  can  be  considered  as  a  thin  region  of  atomic  disarray  where  the 
density  of  atoms  is  slightly  less  than  normal.  Impurity  and  second  phases  are  located  at 
grain  boundaries  [Davidge,  1979].  The  structure  of  the  grain  boundary  in  alumina  was 
studied  by  Carter,  et  al.  (1980)  using  electron  diffraction  and  weak- beam  (or  dark  field) 
imaging  techniques  in  a  Transmission  Electron  Microscope  (TEM).  The  diffraction 
patterns  and  images  show  that  the  structure  of  the  grain  boundaries  is  periodic  and  the 
thickness  of  the  grain  boundaries  is  about  6-9  nm,  which  consists  of  dislocation  arrays  and 
facets. 


4.5  Real-Time  Fracture  Study  of  Ceramics 

In  order  to  perform  the  real-time  fracture  study  of  ceramic  materials,  a  closed-loqp 
control  system  for  the  in-si tu  loading  stage  of  die  Scanning  Electron  Microscope  (SEM)  is 
needed.  For  simplicity,  the  loading  stage  was  designed  outside  the  SEM  chamber  to  test  its 
ability  and  will  be  mounted  inside  the  SEM  chamber  in  the  future. 

Closed-Loop  Control  Loading  Stage 

Fig.  3  (a)  shows  the  tension-compression  (T-C)  stage  which  can  be  mounted  to  the 
SEM  position  stage.  A  single  sliding  brass  block  moves  in  a  machined  aluminum  body 
assembly.  To  implement  rotation  of  the  lead-screw,  a  worm  and  ring  gear(40: 1  ratio)  set¬ 
up  is  utilized.  This  set-up  greatly  reduces  reverse  rotation  and  provides  a  smooth 
application  of  lead-screw  torque.  A  square  shaft  slides  in  a  keyed  slot  running  through  the 
worm  gear,  and  a  90  degree  bevel  gear  system  with  another  square  free-sliding  shaft 
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connects  through  a  Detain™  U-joint  to  a  vacuum  rotational-feed-through  shaft  This  shaft 
can  be  routed  by  the  geared  stepper  motor  for  precise  feed-back  control. 

The  closed  loop  design  must  have  several  properties,  including,  but  not  limited  to, 
accurate  load  or  displacement  increments.  Therefore,  die  objective  is  to  create  an  integrated 
environment  between  all  of  the  hardware  elements  necessary  to  perform  tests.  See  Fig.  3 
(b)  for  a  schematic  representation  of  the  operating  environment. 

To  give  an  overview  of  operating  principles,  knowledge  of  the  involved 
components  and  hardware  is  needed  including: 

1:  Scanning  Electron  Microscope 

2:  Loading  stage  capable  of  tension  and  compression  testing 

3:  Step  Motor  and  Controller  from  American  Precision  Industry 

4:  Galil  Motion  Control  Card 

S:  Crack  Opening  Displacement  (COD)  gage  from  MTS 

6:  1000  lb.  Load  Cell  from  Sensotech 

7:  2-channel  DC  amplifier  from  Ectron 

8.  Microsoft's  Virtual  Basic 

The  Galil  Motion  Control  card  is  the  interface  between  the  computer  program  to 
close  the  control  loop  and  the  step  motor  controller,  which  controls  the  motor,  and  puts 
load  and  displacement  on  die  specimen.  The  card  also  receives  two  transducer  signals  from 
the  loading  stage.  The  COD  gage  returns  the  displacement  between  loading  platforms  and 
the  load  cell  returns  the  approximate  load  on  the  specimen  (see  Fig.  3).  The  actual  closed- 
loop  program  in  the  Galil  card  language  has  been  completed.  Microsoft’s  Virtual  Basic  is 
being  used  to  create  a  driver  program  which  is  more  visual  and  convenient.  A  detail 
discussion  on  the  development  of  tension-compression  loading  stage  is  addressed  in 
Appendix  3. 

Compact  Tension  Specimen  Preparation 

Because  sample  preparation  plays  an  important  role  in  the  ceramic  test  [Rice,  1993], 
the  steps  of  sample  preparation  include  the  following: 

1.  The  alumina  bar  was  cut  into  small  pieces  by  the  diamond  saw,  and  then  the 
Isomet  Low  Speed  Saw  (made  by  Buehler  Ltd.)  was  used  to  section  the  small  piece  into  1- 
4  mm  thickness  samples  30  mm  X  30  mm,  sectioned  by  using  a  low  concentration 
diamond  watering  blade. 

2.  The  Ecomet  Grinder  and  Polisher  was  used  to  grind  and  polish  the  surface.  The 
diamond  root  disc  and  diamond  lapping  film  from  40  micron,  15  micron,  6  micron,  3 
micron  to  0.5  micron  were  used  to  obtain  the  shiny  surface. 

3.  The  holes  on  the  specimen  were  drilled  by  an  Ultrasonic  Disc  Cutter  (made  by 
Gatan  Inc.).  The  disc  cutter  vibrates  a  tubular  tool  at  a  frequency  of  about  26  Kz  against  the 
sample.  Hie  tool  is  immersed  in  a  drop  of  water-based  hard  grit  (usually  SiC  320  grit 
powder)  slurry  placed  on  the  specimen.  This  causes  particles  in  the  slurry  to  inpact  the 
sample  under  the  vibrating  tool  and  erode  away  a  circular  impression. 

4.  Notches  were  sawed  by  an  Isomet  Low  Speed  Diamond  Saw  with  a  thin  diamond 
wafering  blade.  Usually  die  chevron  notch  is  cut 
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5.  After  machining  the  specimens,  the  specimens  were  heated  to  600  °C  for  2  hours 
by  a  Thermolyne  temperature-controlled  furnace  to  release  the  residual  stress  caused  by  the 

machining.  Hie  heating  and  cooling  rate  was  2°C/min. 

6.  Aluminum  was  evaporated  on  the  specimen  surface  by  a  Denton  Vacuum  DV- 
515  to  produce  a  thin  conductive  layer. 

Experimental  Results 

Grain  boundaries  in  ceramic  polycrystals  frequently  constitute  planes  of  reduced 
fracture  resistance;  consequently,  crack  propagation  often  occurs  along  grain  boundaries. 
Although  grain  boundary  or  intergranular  fracture  does  not  necessarily  give  rise  to  material 
weaker  than  materials  that  fail  by  transgranular  fracture,  the  preferential  failure  along  grain 
boundaries  can,  in  some  cases,  enhance  toughness.  To  provide  insight  into  the  toughness¬ 
determining  mechanisms,  both  fracture  paths  and  fracture  surface  were  examined  during 
load  application. 

Fig.  4  and  Fig.  1 1  show  the  area  around  the  notch  tip  before  and  after  the  specimen 
was  broken,  respectively.  Load  was  applied  carefully  by  constantly  detecting  the  notch  tip 
at  high  magnification.  A  crack  started  to  propagate  from  the  notch  tip  (see  Fig.  5)  after  a 
certain  load.  It  can  barely  be  seen  in  low  magnification.  The  profiles  of  the  crack  while 
holding  the  applied  load  are  shown  in  Fig.  6  (a),  (b),  (c),  (d),  (e).  The  micrographs  show 

that  the  crack  started  to  propagate  from  the  notch  tip  to  1 100  pm  in  length.  The  width  of  the 

crack  was  1-2  pm.  It  demonstrated  significant  crack  deflection  with  the  tortuosity  of  the 
crack  path  (Fig.  7  (a)  and  (b)). 

Specific  examples  of  SEM  observations  in  alumina  are  shown  in  Fig.  8  (a)  and  (b). 
Fig.  8  (a)  shows  a  bridging  site  in  the  crack  wake  some  950  pm  behind  the  crack  tip  (also 
shown  in  Fig.  5  (a));  Fig.  8  (b)  shows  another  bridging  site  in  the  crack  wake  some  1000 
pm  behind  the  crack  tip  (also  shown  in  Fig.  5  (a)).  Other  bridging  sites  were  also  observed 
along  the  profile  of  the  crack. 

The  fracture  mode  was  predominantly  intergranular,  and  it  also  was  observed  on 
the  fracture  surface  (Fig.  9).  No  indication  of  a  microcrack-cloud  zone  (also  known  as 
frontal-zone  microcracking)  was  observed  by  SEM  (see  Fig.  10). 

Summary  of  experimental  results: 

1.  Crack-deflection  processes  operate  when  a  crack  interacts  with  microstructurally  related 
features  (e.g.,  weak  interfaces  or  residual  stress  fields)  that  reduce  deviation  from 
planarity.  As  the  crack  deflects  out  of  the  plane  normal  to  the  applied  stress,  the  stress 
intensity  at  the  tip  diminishes,  reducing  the  crack  driving  force  and  improving  the  fracture 
toughness.  K.  T.  Faber  and  Anthony  G.  Evans  found  intergranular  crack-deflection 
toughening  in  silicon  carbide  [Faber,  1983]. 

2.  Frontal-zone  microcracking  was  not  found  around  the  fracture  tip. 

3.  Since  no  energy  was  dispersed  in  front  of  the  crack  tip,  crack  interface  grain  bridging  in 
the  wake  of  the  crack  might  be  another  cause  of  fracture  toughening  [Rodel,  1992]. 

4.  The  failure  mode  was  predominantly  intergranular. 
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4.4  Numerical  Simulation  of  Microfracture 


Based  on  the  microstructural  study  of  ceramic,  it  was  proposed  to  numerically 
model  in  two  dimensions  a  small  specimen  composed  of  a  set  of  grains  and  grain 
boundaries,  as  shown  in  Fig.  1.  Each  grain  can  be  roughly  pictured  as  a  homogeneous, 
anisotropic  material  with  a  multifaceted  surface  consisting  of  several  planes.  When  a 
ceramic  is  loaded  statically  in  tension,  micrographs  show  that  a  dominant  response 
mechanism  is  the  appearance  of  cracks  along  gram  boundaries  while  the  individual  grains 
remain  intact  Therefore  the  mechanical  features  of  the  grain  boundary  itself  must  have  a 
significant  impact  on  the  mechanical  behavior  of  the  ceramic.  It  is  surprising  that  (to  our 
knowledge)  no  attempt  has  been  made  to  characterize  the  mechanical  features  of  the  grain 
boundary  itself  since  it  is  the  grain  boundary  that  appears  to  provide  the  dominant 
characteristics  of  the  mechanical  behavior  of  ceramics.  It  is  generally  recognized  that  the 
nature  of  grain  boundaries  and  microstructure  affect  the  properties  and  behavior  of 
ceramics.  It  is  equally  important  to  understand  ceramic  processing  characteristics  which 
affect  the  nature  of  the  gram  boundaries.  In  general,  the  mechanical  properties  of  ceramics 
depend  on  the  strength  at  grain  boundaries,  intergranular  fracture  due  to  weak  grain 
boundaries  and  transgranular  fracture  due  to  strong  grain  boundaries. 

Grain  boundaries  are  considered  as  thin  regions  of  elastic  damaging  material.  An 
initial  assumption  is  that  each  grain  is  isotropic.  The  grain  boundary  is  supposed  to  be 
restricted  to  mode  I  behavior  in  the  direction  perpendicular  to  the  face  of  the  grain.  This 
response  is  accomplished  by  utilizing  a  continuum  model  with  a  high  shear  modulus  to 
preclude  grain  boundary  sliding,  a  continually  increasing  modulus  for  compressive 
deformation,  and  a  decreasing  stiffness  modulus  as  a  damage  model  to  simulate  the  mode  I 
response  for  a  tensile  stress  normal  to  the  face  of  the  grain. 


Continuum  Damage  Mechanics 

The  nonlinear  behavior  of  brittle  materials  is  more  accurately  represented  as  the 
evolution  of  distributed  microcracks  rather  than  plastic  deformation;  therefore,  continuum 
damage  mechanics  is  the  most  appropriate  approach.  A  phenomenological  model  [Yazdani 
and  Schreyer,  1988]  has  been  developed  based  on  perceived  modes  of  crack  evolution.  The 
value  of  this  work  is  in  its  establishment  of  a  firm  thermodynamical  foundation  for  the 
approach,  and  when  combined  with  plasticity  [Yazdani  and  Schreyer,  1990],  the  model 
provides  excellent  qualitative  and  quantitative  results  for  a  variety  of  load  paths.  Some 
theoretical  implications  of  the  use  of  continuum  damage  mechanics  have  been  explored 
[Schreyer  and  Wang,  1990]  as  a  precursor  to  the  incorporation  of  essential  aspects  of 
microstructural  approaches  [Ju,  1990;  Krajcinovic,  et  al.,  1990]  in  a  sufficiently  simple 
manner  to  retain  the  feasibility  of  performing  numerical  simulations. 

Conventional  damage  consists  of  the  creation  of  voids  and  microcracks  which 
consequently  reduce  the  mechanical  properties  of  a  material,  such  as  bulk  modulus  and 
shear  modulus.  Although  damage  is  an  anisotropic  phenomenon,  the  initial  approach  will 
be  to  assume  isotropy  and  consider  the  pan  of  the  damage  reflected  through  the  bulk 
modulus.  A  preliminary  theory  consistent  with  thermodynamics  suggests  that  volumetric 
strain  is  a  suitable  measure  of  damage. 

If  the  thermal  effects  are  ignored,  the  internal  energy,  U,  is  assumed  to  be  a 
function  of  the  total  strain  tensor,  e,  the  "internal"  variable  consisting  of  the  permanent 

strain  tensor,  eP,  and  the  fourth-order  elasticity  tensor,  E.  With  the  assumption  of  linear 
elasticity,  U  is  taken  to  be 
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«  U  *  j  (e-«P):E:(e-«P) 

the  application  of  the  first  and  second  (Clausius-Duhem  inequality)  laws  of 
thermodynamics  leads  to  the  constitutive  relation  for  stress,  a, 

#  a  «  •  ■  E:(e-eP) 

and  the  dissipation  inequality 

T  *  T  •  •  |  • 

#  0  °r  attf-  j(e-«P) :  E:  (e-eP)^O 

Suppose  the  damage  process  is  parameterized  through  die  use  of  a  parameter,  o>,  which  is 

monotonically  increasing  (<o£  0).  The  evolution  equations  for  damage  can  be  given  as 
A  follows: 


E  « -  o>R(E,  e,  eP  )  e^ta  m(E,  e,  eP) 


in  which  the  dependence  of  the  response  functions,  R  and  m,  is  shown.  Since  0,  the 
dissipation  inequality  becomes 

D=  j  (e-eP):R:(e-«P)  +  s:m  2: 0 


which  is  satisfied  if  R  is  positive  definite  (conventional  definition  of  damage)  and  if  m 
forms  an  acute  angle  with  s. 

Suppose  a  damage  function,  f,  is  defined  such  that  damage  occurs  when  the  state  falls  on 
the  damage  surface,  f  ■  0,  and  for  which  the  dissipation  inequality  is  satisfied 
automatically.  The  inequality  is  met  if 

f  *  D  -  g2(E,  e,  eP) 

Then  when  f  <  0  damage  is  not  occurring,  and  f  >  0  is  not  a  physical  state. 

Let  I  and  i  denote  the  fourth  order  and  second  order  identity  tensors,  respectively.  Define 
the  spherical  and  deviatonc  projections  to  be 

P*P-ji®i  pd.I-psp 

Then  P*P :  P*P  -  PSP  pd  :  pd  «  pd  PSP  :  pd  ■  0 

ad  *  pd;S  ed  =  pd:e 
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*SP  «  psp:s  .  -Pi  e*p  «  P*P:e  -  5  Cy  I 

in  which  P  *  -{i:s)/3  is  the  mean  pressure,  ey-Cytl  is  the  volumetric  strain,  and  sd  and  ed 
denote  the  stress  and  strain  deviators,  respectively.  With  the  use  of  the  projection 

operators,  the  isotropic  elasticity  tensor  is  E  *  3KPsP+2GPd  in  which  K  is  the  bulk 
modulus  and  G  is  the  shear  modulus.  Both  parameters  may  change  with  damage. 

For  the  initial  numerical  investigation  performed  for  this  study  we  chose  g  to  be  a 
constant,  and  R  to  be  the  initial  isotropic  elastic  tensor  Eg,  and  m  to  be  the  spherical  part  of 

the  stress  tensor  (tensile).  The  result  is  a  simple  scalar  isotropic  damage  model  in  which 
the  bulk  and  shear  moduli  decay  simultaneously  with  co  which  is  directly  related  to  the 
volumetric  strain  (positive).  A  slightly  more  general  model  was  then  introduced  in  which 
R  is  an  isotropic  tensor  not  proportional  to  E0.  The  result  is  a  two  parameter  isotropic 

damage  model  in  which  the  bulk  and  shear  moduli  deteriorate  at  different  rates  with  an 
increase  in  volumetric  strain. 


Damage  Model 

For  the  3-D  problem,  the  material  stiffness  matrix  [E]  is  6  by  6.  In  the  most 
general  case  of  anisotropy,  [E]  contains  21  independent  coefficients.  AM  coefficients 

change  as  the  damage  parameter  o)  changes.  That  is  why  continuum  damage  mechanics  is 
so  complicated  even  for  the  2-D  problem.  In  the  engineering  sense,  many  situations  can 
be  simplified  or  idealized.  For  an  orthotropic  material,  the  material  displays  extreme  values 
of  stiffness  in  mutually  perpendicular  directions,  and  [E]  contains  9  independent 
coefficients;  for  an  isotropic  material,  the  material  behavior  is  the  same  in  all  directions, 

and  material  properties  are  commonly  expressed  as  two  independent  coefficients,  E  and  p , 
or  K  and  G. 


In  the  2-D  problem,  two  different  situations  should  be  considered;  the  plane  strain 
and  the  plane  stress  problem.  For  the  plane  strain  isotropic  problem,  [E]  can  be  expressed 
as 


[E] 


(l+p)(l-2p) 


"  1-11 
11 

11  n 

3(l-p)K 

m  0] 

p  u 
1-p  0 

n  l-2p 

m 

G+li) 

3pK 

(l+li) 

3(1-P)K  Q 

0 

0  — - 
2 

(1+lD 

(l+li) 

4m 

0 

0  G- 

and  die  volumetric  strain  is 

evm£xx  +  eyy 
For  plane  stress. 
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and  the  volumetric  strain  is 

ev*exx+eyy+eBt 

*^(8xx+eyy) 

where  E  and  |i  are  Young's  modulus  and  Poisson's  ratio,  respectively,  and  K  and  G  are 
bulk  modulus  and  the  shear  modulus,  respectively. 

It  is  assumed  that  the  damaging  regions  are  confined  to  the  grain  boundaries,  and 
the  grains  are  elastic.  For  simplicity,  the  grains  are  assumed  to  be  isotropic  elastic,  and 
the  proportional  scalar  isotropic  damage  was  used  in  the  grain  boundaries,  i.e.,  the  bulk 

Kgb  and  shear  modulus  Ggb  of  the  grain  boundaries  decreased  when  their  volumetric 
strains  reach  a  certain  volume.  The  superscript  GB  denotes  the  grain  boundary.  A 

damage  parameter  oo  is  introduced: 

Kgb«(1-©)k£B 

Ggb  *  ( 1  -  cd  )  G^b 

If  cd  *  0,  no  damage  has  occurred;  if  to  =  1,  complete  damage  has  occurred.  K^B  and 
Gf  are  initial  values  of  K°B  and  GGB,  respectively.  KGB  (and  GGB )  reduces  from 
Kf  (and  G0B )  toOasco  goes  from  0  to  1. 

The  evolution  of  the  damage  parameter  co  under  tension  is  proposed  to  be  (Fig.  12) 


od«  0>(eo,ev,k) 
0 


< 


cv  '  co* 


1-  exp  ( -  (— g- — )  k)  ev  *  ec 


where  ey  is  volumetric  strain,  and  Cq  and  k  are  parameters  which  can  be  determined  from 
data 

•  The  Young's  modulus  E  GB  is 
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Dynamic  Relaxation 

Nonlinear  structural  analyses  include  materially  nonlinear  problems  such  as 
nonlinear  constitutive  equations  with  small  deformations,  geometrically  nonlinear 
problems  normally  associated  with  buckling  or  a  combination  of  both  types  of 
nonlinearities  [Chen,  and  Schreyer,  1990].  The  nonlinear  analysis  in  the  finite  dement 
method  can  be  expressed  as 

[K]  |U)  *  {  F  | 

where  [K]  is  the  structure  stiffness  matrix, 

{11}  is  the  vector  of  nodal  degree  of  freedom,  and 
(F)  is  the  vector  of  nodal  farads. 

[K]  and  {F]  are  regarded  as  dependent  on  {U}.  The  schemes  for  nonlinea  problems  are 
based  on  step-by-step  load  incrementation  and  an  iteration  procedure  to  correct  the 
linearization. 

In  our  finite  element  code,  dynamic  relaxation  (DR)  is  used  as  an  explicit 
incremental  method  for  solving  static  problems.  The  DR  method  is  based  on  the  fact  that 
the  static  solution  is  die  steady  state  pan  of  die  transient  response  for  a  temporal-step  load. 
This  method  is  especially  attractive  for  problems  with  highly  nonlinear  geometric  and 
material  behaviors,  which  includes  limit  points  and  regions  of  very  soft  stiffness.  Use  of 
die  dynamic  relaxation  method  makes  the  software  simple.  Since  the  method  is  explicit,  it  is 
unnecessary  to  form  a  global  stiffness  matrix;  therefore,  much  less  computer  storage  is 
needed  than  with  implicit  procedures  such  as  Newton-type  methods  [Underwood,  1983, 
Gentle,  and  Xie,  1992].  See  Appendix  1  for  the  DR  algorithm. 

CoMeT 

The  whole  modeling  process  was  accomplished  by  coupling  the  nonlinear  dynamic 
relaxation  analysis  program  with  the  program  CoMeT  (Computational  Mechanics  Toolkit), 
and  an  interactive  graphical  shell  for  integrated  computational  mechanics.  For  each 
iteration,  the  new  damage  values  of  K  and  G  were  calculated  in  terms  of  the  current  volume 
strain  in  each  Gaussian  point  and  compared  with  the  former  damage.  After  one  step  was 
finished,  a  data-base  file  was  output  that  could  be  processed  by  CoMeT. 


Mesh 


Based  on  crystal  structure  and  micrographs  of  the  alumina,  the  shape  of  each  grain 
was  hypothesized  as  hexagonal  and  modeled  with  six  three-node  triangular  elements;  the 
uniform  grain  boundaries  are  distributed  between  the  grains,  which  was  represent  with  a 
small  three-node  triangular  element  and  a  four-node  quadratic  element  (  Fig.  1).  For 
simplicity,  a  small  mesh  with  a  grain  size  of  10  mm  was  generated,  and  the  thickness  of 

•  the  grain  boundaries  was  0.5  mm  (Fig.  1).  The  aspect  ratio  erf  the  grain  boundary  element 
was  10:1.  An  element  performs  best  if  its  shape  is  compact  and  regular.  According  to 
Code  [1989],  an  element  tends  to  stiffen  and  lose  accuracy  as  its  aspect  ratio  increases, 
although  specific  details  are  not  given.  Here,  we  show  that  good  results  are  achieved  even 
with  a  large  aspect  ratio. 

•  Numerical  Results 

To  illustrate  the  proposed  damage  model,  two  sample  problems  were  considered; 
one  was  a  small  problem  with  19  elements,  another  was  a  large  problem  with  161 
elements.  The  model  problems  consist  of  a  specimen  in  which  the  displacement  on  one 
boundary  is  applied  uniformly  in  the  y  direction  under  the  assumption  of  plane  stress.  One 
9  boundary  was  traction  free  and  the  other  two  boundaries  had  fixed  displacement  in  the  x 

and  y  directions,  respectively.  In  the  grain  regime,  the  material  is  assumed  to  be  isotropic 

as  defined  through  Young's  modulus,  E®,  and  Poisson's  ratio,  pG.  Material  parameters, 
which  are  considered  to  be  representative  of  alumina  ceramic,  were  chosen  as  follows:  E® 
*  372  GPa,  and  |iG  =  0.22.  Equivalently,  bulk  modulus  KG=  228  GPa,  shear  modulus 

•  Gc  =  152  GPa,  were  used  for  the  grain.  The  initial  material  properties  of  grain  boundaries 
affect  the  computational  results. 

Small  Mesh  Model  Problem 

The  finite  element  mesh  and  the  boundary  conditions  used  in  the  analysis  are  shown 

•  in  Fig.  13.  First,  the  initial  material  properties  of  grain  boundaries  used  in  the  numerical 

OB  G 

analysis  were  assumed  as  the  same  values  as  that  for  the  grains,  i.e.,  K  0  *  K  *  228 
GPa,  GGB  =  G°  =  152  GPa;  equivalently,  EGB  =  372  GPa,  and  p  GB  «  0.22.  The 

undeformed  mesh  is  show,  in  Fig.  13.  The  deformed  meshes  at  the  20^  step  and  30^ 

•  step  are  shown  in  Fig.  14(M)  and  14(b  ).  The  larger  deformation  that  occurred  in  the 
shadow  elements  showed  that  the  microcrack  started  to  nucleate  and  propagate  along  the 
grain  boundaries. 

Fig.  14  (a)  shows  the  crack  nucleating  at  the  central  grain  boundary  (marked), 
and  crack  propagation  along  the  grain  boundaries  is  shown  in  Fig.  14  (b).  The 

•  displacement  and  the  total  corresponding  nodal  force  curve  is  illustrated  in  Fig.  14  (c). 
The  nonlinear  behaviors  (hardening  and  softening)  of  ceramics  under  uniaxial  tension  can 
be  seen.  The  filled  circles  indicate  the  loading  steps. 


Large  Mesh  Sample  Problem 

Hie  large  mesh  sample  problem  includes  12  grains  and  the  number  of  elements  is 
161.  The  mesh  is  shown  in  Fig.  15  (a).  The  features  of  the  proposed  model  are  shown 
below: 
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1.  Microcnck  nucleation: 


The  undefocmed  mesh  with  161  elements  is  shown  in  Hg.  IS  (a).  With  E"  * 

E°,  the  deformed  mesh  at  the  23IX*  step  is  shown  in  Fig.  15(b).  It  can  be  seen  (Fig.  IS 
(b))  that  die  grain  boundaries  (in  shadow)  show  a  larger  deformation  than  the  other 
elements,  i.e.,  die  microcncks  have  been  created  on  the  grain  boundaries.  Hg.  IS  (c)  gives 

the  displacement  and  corresponding  total  nodal  force  curve  for  die  condition  of  ■  E°. 

2.  Dick  propagation: 

As  the  displacement  is  continuously  applied,  the  deformed  mesh  shows  that  the 
tmcrocracks  connect  together  and  propagate.  Hg.  16  (a),  (b),  (c)  and  (d)  show  results  at 

the  25th,  26th,  27*  and  291*1  load  steps.  From  the  deformed  mesh,  crack  nucleation  and 
propagation  along  the  grain  boundary  can  be  seen  dearly. 

3.  Microcrack  and  void  effects: 

The  program  could  also  simulate  dements  as  tmcrocracks  and  voids  when  the  initial 

mechanical  properties  of  grain  boundaries,  like  E^B,  were  given  as  very  small.  One 

dement  was  simulated  as  a  microcrack  as  shown  in  Hg.  17  (a)  (marked  dements).  The 
deformed  mesh  is  shown  in  Hg.  17  (b)  and  (c).  Compared  with  Hg.  IS  (b),  Hg.  17 
(b)  shows  that  the  crack  starts  to  propagate  from  the  weakened  element  (marked  element), 
i.e.,  the  crack  propagates  from  the  microcracks.  Hg.  17  (c)  shows  die  difference 
between  the  force  vs.  displacement  curves  with  and  without  existing  microcracks.  The 
existing  microcrack  reduced  the  strength  of  the  ceramic  and  exhibited  more  brittle  material 
properties  than  the  ceramic  without  nucrocracks. 

4.  Inclusion  effect: 

Fig.  16  (d)  shows  that  the  main  crack  propagates  along  die  grain  boundaries  in  the 
upper  part  of  the  specimen.  The  inclusion  was  introduced  in  die  marked  element,  i.e.,  a 
larger  value  of  stiffness  was  used  in  die  marked  element  (H|.  18(a)).  Horn  die  deformed 
mesh,  the  position  of  the  main  crack  was  changed  due  to  die  introduced  inclusion  element, 
as  can  be  seen  in  Hg.  18  (b)  and  (c).  Compared  with  Rg.  16,  die  inclusion  deterred  the 
crack  nucleation  ana  propagation  which  would  ocherwwue  run  through  it  before.  Hg.  18 
(d)  shows  die  difference  between  die  force  vs.  displacement  curves  with  and  without  an 
inclusion.  The  strength  of  the  ceramic  was  enhanced  by  introducing  an  inclusion,  but 
ductility  was  decreased. 

5.  Effects  of  mechanical  properties  of  the  grain  boundary  on  macro-behavior 

In  the  above  computations,  the  initial  mechanical  property  of  die  grain  boundary 
was  assumed  to  be  the  same  as  thegrain.  Generally,  the  grain  boundary  is  weaker  than  the 
grain  for  intergranular  fracture.  Hg.  19  gives  the  displacement  and  corresponding  total 

nodal  force  curve  for  the  case  of  -  E°  and  E®8  «  0.75  E°.  Mechanical  properties  of 

die  grain  boundary  have  a  significant  effect  on  die  macroscopic  behaviors.  The  weaker  the 
grain  boundaries,  die  less  the  strength  and  ductility  of  ceramic.  The  failure  mode  is  the 
same  as  that  shown  in  Hg.  16. 
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Voronoi  diagram 

One  of  the  major  obstacles  in  numerical  modeling  of  ceramic  is  die  geometrical 
complexity  in  grain  microstructure.  Since  many  properties  of  ceramics  are  determined  by 
their  grain  microstructure,  prediction  and  control  of  microstructural  characteristics  are 
important  So  far,  a  number  of  computer  simulations  have  been  made  for  microstructures 
[Anderson,  et  al,  1984, 1989,  Wejchert,  1986]. 

One  of  the  methods  for  simulating  the  grain  geometry  and  topography  is  die  model 
known  as  the  Voronoi  diagram.  Given  a  number  of  points  in  die  plane,  the  Voronoi 
diagram  divides  the  plane  according  to  the  nearest-neighbor  rule:  Each  point  is  associated 
with  the  region  of  the  plane  closest  to  it;  (Fig.  20  (a»  [Auxenhammer,  1991]. 

A  program  for  generating  the  Voronoi  diagram  written  in  C  language  is  available  in 

.1  •  •  «  «w  _  _ _ j  *• _ « _ j  a _ am* _ ^ _ l*  / _ *  •  *L_ 


of  the  vertices  erf  polygons  and  line  connections  was  transferred  to  die  Sun  Station,  then  the 
Voronoi  diagram  was  redisplayed  by  CoMeT  as  shown  in  Fig.  20  (b). 

For  finite  element  modeling,  a  thin  layer  between  grains  was  generated  as  a  grain 
boundary.  The  grain  was  meshed  with  3-node  triangular  elements,  and  the  grain  boundary 
was  meshed  with  4-node  quadratic  elements  and  3-node  triangular  elements  (Fig.  21). 

Scalar  isotropic  damage  mechanics 

For  simplicity,  isotropic  damage  mechanics  was  used  first.  As  we  mentioned 
before,  the  microcracks  started  to  nucleate  at  the  grain  boundaries,  and  then  propagated 
along  the  grain  boundary,  i.e.,  an  intergranular  fracture  as  shown  in  Fig.  22  (a),  (b),  (c), 
and  (d).  The  grain's  rotation  can  be  observed  in  the  deformed  mesh. 

Two-parameter  isotropic  damage  mechanics 

The  grain  boundary  is  restricted  to  mode  I  behavior,  perpendicular  to  the  face  of  die 
grain.  This  response  is  accomplished  by  utilizing  a  continuum  model  with  a  higher  shear 
modulus  to  preclude  grain-boundary  sliding,  a  continually  increasing  modulus  for 
compressive  deformation,  and  a  decreasing  stiffness  modulus  as  a  damage  model  to 
simulate  the  mode  I  response  for  a  tensile  stress  normal  to  the  face  of  the  grain. 

The  material  stiffness  matrix  [E]  can  be  expressed  as 


3(1 -p)K  3nK  Q 

(1+ji)  (1+p) 

3uK  3(l-ii)K  Q 
(1+JA)  (l4|i) 

0  0  G 


(for  plane  strain) 
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(for  plane  stress) 


where  K  and  O  are  the  bulk  modulus  and  shear  modulus,  respectively  in  general  cases.  The 
material  stiffness  matrix  [E]  changes  as  K  and  O  change,  and  Poisson's  ratio  p  remains 
constant  for  alumina:  p  -0l22  -  0.23. 


Two  damage  parameters  00]  and  cifr  were  introduced: 

K°b.(1-<d1)IcSB 

G®.(1-«12)0“ 


and  0)]  >0)2  (see  Fig.  23). 


Using  the  two-parameter  isotropic  damage  model,  the  computational  results  are 
shown  in  Kg.  24  (a),  (b),  (c),  and  (d).  The  microcracks  started  to  nucleate  at  the  grain 
boundaries,  and  propagated  along  the  grain  boundary  until  finally  a  macrocrack  was 
formed. 

Kg.  23  shows  die  difference  of  the  load  and  displacement  curves  between  the  two 
types  of  damage  modes.  The  grain  boundaries  with  a  higher  shear  modulus  prevented  the 
rotation  of  grams,  but  the  strength  and  ductility  of  ceramics  were  reduced. 


Summary  of  numerical  investigation: 


1.  Based  on  the  continuum  daman  mechanical  model,  microcrack  nucleation,  crack 
propagation,  microcrack  and  void  effect,  and  inclusion  effect  were  simulated  numerically 
by  a  small  specimen  composed  of  a  set  of  grains  and  grain  boundaries  in  two  dimensions. 

2.  The  grain  size  can  be  randomly  generated  by  the  Voronoi  Diagram,  so  that  the 
computational  mesh  can  be  more  realistic.  The  tortousity  of  the  crack  path  and  crack 

•  bridging  were  observed  by  the  numerical  results. 

3.  Because  of  die  computational  limitation  of  die  work  station,  die  size  of  mesh  cannot  be 
generated  large  enough  to  compare  the  actual  size  of  specimens,  but  many  fracture 
mechanisms  can  still  be  observed  m  numerical  investigations. 


4.5  Summary  and  Conclusion 
Experimental 

It  is  generally  recognized  that  the  nature  of  grain  boundaries  and  nricrostructures 
affect  the  properties  and  behavior  of  ceramics.  It  is  important  to  understand  ceramic 
processing  characteristics  which  affect  the  nature  of  the  grain  boundaries.  Different 
#  processing  techniques  have  been  explored  to  produce  the  desired  microstnictures  for 

ceramics.  Several  processing  techniques  have  been  explored  to  produced  the  desired 
microstructure  alumina.  It  was  found  that  Two-Stage  sintering  and  MgO-doped  AljOj 
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can  be  combined  to  produce  a  high  performance  alumina.  The  result  is  better  alumina  with 
I  an  improved,  more  uniform  and  homogenous  microstructure. 

Examination  of  micros  tructures  has  indicated  that  mechanical  properties  depend  on 
the  population  and  location  of  pores  and  grain  size  distribution.  The  fracture  resistance  of 
ceramics  can  be  affected  by  micros  true  tural  variables  such  as  porosity  and  grain  size 
distribution.  Larger  grain-size  ceramics  are  usually  more  prone  to  microfracture.  Porosity 
1  located  as  grain  boundaries  degrades  the  fracture  resistance  in  proportion  to  die  volume 

concentration  of  voids;  however,  voids  within  grains  can  enhance  toughness  during 
tnuisgranular  fracture. 

An  image  analyzer  and  an  image  editor  (NIH  IMAGE  1.41)  were  used  for 
quantifying  the  grain  size  distribution  from  micrographs.  The  micrographs  were  redrawn 
g  by  using  a  transparent  material  so  that  the  grains  and  grain  boundaries  could  be  clearly 

identified.  The  image  was  obtained  by  scanning  the  redrawn  micrographs  into  a  computer. 
A  digital  readout  displayed  the  area  for  each  grain  from  the  images.  Using  a  ruler,  the 
multiplying  factor  was  obtained  by  converting  the  picture-point  units  to  units  of  square 

micrometers.  From  this  area,  the  diameterOim)  of  the  equiarea  circle  was  computed.  The 
^  diameter  is  used  as  the  grain  size  measure. 

In  each  micrograph,  different  magnifications  were  used.  Using  the  statistical 
program  as  developed  through  this  research,  we  obtained  the  final  result  of  the  frequency 
grain  size  distribuuon  for  each  specimen  using  different  sintering  processes.  The  grain  size 
distributions  are  given  in  Fig.  26.  The  numbers  of  grains  for  each  grain  size  distribution 
were  100,  250, 124, 166,  222,  and  215  corresponding  to  sintering  processes  #1  through 
i  #6,  respectively.  The  frequency  grain  size  distribution  shifts  toward  the  coarse  end  of  the 

size  spectrum  as  a  function  of  sintering  time  and  temperature.  Similar  conclusions  can  also 
be  made  for  void  size  distribution  where  average  void  size  shifted  to  larger  void  size  as  a 
function  of  sintering  time  and  temperature  as  shown  in  Figure  27.  It  is  known  that,  during 
sintering,  very  high  temperatures  and  long  periods  of  time  lead  to  exaggerated  grain 
^  growth.  The  grain  size  distribution  can  be  approximately  fitted  to  a  Lognormal  Distribution. 

The  X-ray  diffraction  method  was  used  to  determine  the  shape  and  size  of 
alumina's  unit  cell.  The  data  demonstrate  that  the  alumina  has  a  hexagonal  crystalline 
structure.  In  addition  to  the  X-ray  diffraction  scan,  the  Back-Reflection  Laue  method  was 
used  in  qualitative  texture  analysis  to  determine  the  distribution  of  grain  orientation  in  the 

sintered  alumina.  Crystal  sizes  of  more  than  10  pm  were  found  to  be  present  and  the 

*  alumina  crystals  were  randomly  oriented. 

One  of  the  objectives  of  this  research  is  to  study  the  fracture  process  of  the  ceramic 
in  real  time  and  to  measure  die  strain  field  in  the  vicinity  of  the  crack  including  die  bridging 
zones.  A  fast-scanning  electron  microscope(FSEM)  for  dynamic  microscopy  applications 
f  was  used  to  capture  the  fracture  events  in  die  ceramic.  This  equipment  is  particularly 

*  suitable  for  microstructural  studies  for  dynamically  loaded  materials,  as  it  captures  images 
at  a  high  speed.  The  SEM  chamber  was  also  modified  to  accommodate  an  in-situ  tension- 
compression  loading  device  to  fracture  ceramics.  Several  problems  have  arise  including  die 
inadequacy  Of  force  to  achieve  ceramic  fracturing  and  excessive  noise  that  was  produced  by 
step  motor  instrumentation.  Many  improvements  have  been  made  such  as  application  of  an 

_  adequate  force  driven  by  a  microstepper  motor  and  closed-loop  control  capability  by  using 

*  filtering  techniques. 

Numerical 


i 
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Continuous  damage  mechanics  is  an  appropriate  approach  for  analyzing  the 
nonlinear  behavior  of  qu&si-brittle  materials.  It  was  proposed  to  numerically  model  in  two 
dimensions  a  small  region  composed  of  a  set  of  gnuns  and  grain  boundaries.  Each  grain 
was  modeled  as  an  elastic  material  while  the  grain  boundary  was  considered  to  be  a  thin 
region  of  elastic  damaging  material.  The  randomly  distributed  grain  geometry  was 
simulated  by  a  computer  program  which  produced  the  Voronoi  diagram.  Both  scalar  and 
two-parameter  isotropic  damage  models  were  used  in  the  computation. 

Microcrack  nucleation,  crack  propagation,  microcrack  and  void  effect,  and 
inclusion  effect  were  simulated  numerically  on  a  small  specimen  composed  of  a  set  of 
grains  and  grain  boundaries  in  two  dimensions.  The  tortuosity  of  the  crack  path  and  crack 
bridging  were  observed  through  the  numerical  results. 

Numerical  examples  have  also  demonstrated  the  following  observations: 
Microcracks  always  started  to  nucleate  cm  the  grain  boundary  first  if  the  grain  boundary 
was  weaker  than  the  grain.  If  a  microcrack  already  existed,  then  the  crack  started  to 
propagate  from  the  microcrack  along  the  grain  boundary.  The  existing  microcracks 
weakened  the  ceramic  and  exhibited  a  more  brittle  failure  mode  as  shown  in  Fig.  17(d). 
When  an  inclusion  existed  in  the  grain  boundary,  it  deterred  the  crack  propagation  through 
it,  and  the  crack  initiated  from  the  other  weaker  site.  Inclusions  enhanced  the  strength  of 
ceramic,  but  reduced  its  ductility  (Fig.  18(d)).  Mechanical  properties  of  the  grain 
boundary  have  a  great  effect  on  macroscopic  behaviors.  The  strength  and  ductility  of 
brittle  ceramic  decreased  as  the  strength  of  dm  grain  boundary  decreased. 


5.  Proposed  Extension 

Experimental 

Most  traditional  ceramics  are  monolithic  (single  phase),  such  as  alumina,  and  have 
very  low  fracture  toughness.  Because  they  do  not  yield  to  loading  pressure,  monolithic 
ceramics  behave  as  ideally  brittle  materials,  and  a  propagating  crack  needs  only  to 
overcome  the  surface  energy  of  the  material.  The  new  generation  of  ceramics,  however, 
includes  multiphase  materials  and  ceramic  composites  that  have  vastly  improved  toughness. 
The  micromechanisms  that  lead  to  improved  fracture  resistance  in  modem  ceramics  include 
microcrack  toughening,  transformation  toughening,  ductile  phase  toughening,  fiber 
toughening,  and  whisker  toughening.  Modem  ceramics  also  display  much  larger  ductility 
than  monolithic  ceramics. 

Ceramic  composites  are  considered  an  enabling  technology  for  advanced  aircraft 
and  space  propulsion  engines  and  space  power  systems.  Catastrophic  fracture  remains  an 
issue.  Our  efforts  are  toward  the  in-depth  understanding  of  the  toughening  mechanism 
under  bod)  static  and  dynamic  loading  of  ceramic  composites. 

In  particular,  our  unique  experimental  capabilities  include  a  Fast  Scanning  Electron 
Microscope  (FSEM),  for  which  Dr.  Wang,  the  co- principal  investigator,  received  a  patent 
in  October,  1993.  Specification  and  capabilities  of  the  FSEM  is  described  in  Appendix  4. 
Through  this  research  we  have  also  developed  an  in-si tu  dynamic  closed-loop  controlled 
tension -compression  loading  device.  This  allows  (me  to  observe  in  real-time  the  fracture 
and  toughening  mechanisms  of  reinforced  ceramic-matrix  composites  materials  subjected 
to  static,  cyclic  dynamic  loading. 

Future  research  will  include  the  study  of  toughening  mechanisms  of  compact 
tension  specimens  of  alumina  and  SiC-whisker-reinforced  alumina  using  FSEM  and  an  in- 
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situ  loading  device  inside  the  FSEM  chamber.  Measurements  will  include  load,  fracture 
)  opening  displacement,  crack  extension,  crack-growth  rate  for  cyclic  loading,  number  of 

cycles  to  failure  and  stress  intensity  factor.  Microstructure  observation  will  include  real¬ 
time  fracture  propagation  micrographs  in  a  video  set-up.  From  the  micrographs  one 
should  be  able  to  distinguish  the  fracture  mode  (intergranular  or  transgranular)  and 
toughening  mechanisms  due  to  static  and  cyclic  fatigue  loading. 

I  Theoretical 

A  detailed  understanding  of  ceramics  requires  constitutive  models  for  grains  and 
grain  boundaries,  and  numerical  simulations  of  the  failure  of  grain  boundaries;  the 
resultant  evolution  of  microcracks  can  be  seen  as  crack  branching  and  the  ultimate 
development  of  a  macrocrack.  This  research  has  provided  a  good  initial  analysis  which 
I  shows  that  a  combination  of  continuum  damage  as  a  means  for  simulating  material  failure, 

the  use  of  Voronoi  diagrams  to  construct  realistic  grain  topologies,  and  a  robust  solution 
algorithm  for  computations  beyond  the  peak  force  can  begin  to  show  essential  aspects  of 
the  failure  process. 

Verification  of  the  constitutive  equation  used  for  the  grain  boundary,  and  even  the 
|  grains,  remains  a  difficult  problem.  Since  the  behavior  is  highly  inhomogenous,  and 

because  of  the  small  size  of  a  typical  grain,  it  is  hard  to  imagine  that  direct  correlation  with 
experimental  behavior  will  ever  be  possible.  A  much  more  realistic  approach  is  to 
postulate  reasonable  constitutive  equations  for  grains  and  grain  boundaries,  and  then 
determine  the  response  of  a  specimen  subjected  to  boundary  conditions  analogous  to  those 
applied  experimentally.  Comparisons  of  predicted  load-displacement  curves  and  failure 
g  modes  with  experimental  data  can  then  provide  a  means  for  an  indirect  verification  of 

constitutive  equations. 

If  quantitative  predictions  are  to  made,  a  significant  enhancement  to  the 
computational  approach  is  required.  For  example,  three-dimensional  simulations  must  be 
performed  on  specimens  containing  a  significantly  larger  number  of  grains  if  the  effects  of 
artificial  constraints  imposed  by  boundary  conditions  is  not  to  dominate  the  solution.  Even 
g  with  powerful  computers,  three-dimensional  simulations  must  be  limited  in  scope  if 

solutions  are  to  be  t&tained  at  reasonable  cost.  Every  effort  must  be  made  to  reduce  the 
complexity  of  the  problem  if  such  simulations  are  to  be  feasible. 

Since  the  cunent  finite  element  study  involving  grain  boundaries  in  ceramics  is 
meant  to  be  rally  a  first  step  to  show  that  the  approach  is  plausible,  no  attempt  was  made  to 
incorporate  features  which  might  be  required  for  larger  two-dimensional  and  threc- 
g  dimensional  simulations.  An  important  procedure  would  be  to  constrain  out  modes  that  are 

not  essential  for  the  problem  at  hand,  but  that  can  cause  ill  conditioning  or  excessively 
small  steps  for  dynamic  relaxation.  For  example,  it  is  highly  unlikely  that  tending  modes  in 
the  four-node  quadrilateral  element  are  important  for  the  grain  boundary.  Furthermore,  one 
might  want  to  study  the  extreme  case  of  rigid  grains  (modeled  through  the  use  of 
constraints)  so  that  the  pain  boundary  response  associated  with  a  larger  number  of  grains 
g  could  be  analyzed,  especially  in  three  dimensions. 

In  light  of  its  importance,  it  is  surprising  that  efficient  algorithms  for  the  routine 
incorporation  of  constraints  in  linear  algebraic  equations  is  not  a  subject  widely  studied.  As 
part  of  this  project  a  step  was  taken  to  show  mat  multiple  constraints  could  be  invoked 
directly  without  a  major  modification  to  the  governing  matrix,  and  both  static  and  dynamic 
examples  were  given  for  one-dimensional  problems  [Schreyer  an  Parsons,  1994].  The 
I  results  are  given  in  Appendix  5. 

It  is  proposed  that  this  approach  be  extended  to  two  and  three  dimensions  for 
application  to  the  study  of  grains  and  grain  boundaries.  Another  natural  application  would 
be  to  the  study  of  fibers  in  ceramics  where  the  fibers  could  be  initially  modeled  as 
inextensible  constraints.  Fiber  breakage  could  then  be  simulated  as  a  release  of  the 
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constraint.  There  are  other  potential  uses  for  such  an  algorithm  in  connection  with 
constitutive  equations  in  which  constraints  of  zero  values  for  certain  components  of  stress 
after  material  failure  could  be  imposed 

In  summary,  rapid  advances  in  the  use  of  a  fast-scanning  electron  microscope  are 
providing  unique  experimental  data.  If  these  results  are  to  be  utilized  to  their  greatest  extent, 
numerical  simulations  are  necessary  to  provide  infonnatku  at  the  level  of  grains  and  grain 
boundaries.  Three-dimensional  simulations  are  the  only  realistic  approach,  but  this  will 
require  innovative  numerical  procedures  to  obtain  solutions.  The  efficient  use  of  constraints 
in  conjunction  with  constitutive  equations  based  on  continuum  damage  mechanics  will  be 
an  important  aspect  of  such  a  study. 
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Figure  1.  The  idealized  mesh  of  die  alumina;  die  zooming  of  a  grain  (6  triangle  elements) 
and  grain  boundary  (6  triangle  elements  and  6  quadratic  elements). 
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S  —  Specimen 
G  •—  Worm  Gear 
P  — •  Loading  pins 
C  — -  Load  cell 
LP  —  Loading  platforms 
B  —  Base  of  loading  stage 


Figure  3  (a)  Schematic  of  die  loading  stage 


Figure  3  (b)  Schematic  illustration  of  close-loop  control  of  loading  device. 
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Figure  6  (a)  The  profile  of  the  crack  (crack  propagated  from  upper  part). 
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Figure  6  (c)  The  profile  of  the  crack  (crack  propagated  from  upper  part). 
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Figure  7  (a)  The  tortuousity  of  the  crack  in  the  wake  of  the  crack. 


Figure  7  (b)  The  tortuousity  of  die  crack  in  die  wake  of  the  crack. 
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Figure  8  (b)  The  grain-bridging  site  1000  mm  behind  the  crack  tip. 


Figure  9  Micrograph  of  the  fracture  surface  of  Alumina 


Figure  10  The  crack  tip  at  higher  magnification. 
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Figure  12  (a)  The  proposed  evolution  of  damage  parameter  (O  under  tension. 


Volumetric  Strain 


Figure  12  (b)  The  relative  bulk  and  shear  moduli  decrease  as  volumetric  strain  increases 
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Figure  13  The  finite  element  mesh  and  the  boundary  conditions 
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Figure  14  (a)  The  deformed  mesh  at  die  20th  loading  step 


Figure  14  (b)  The  deformed  mesh  at  the  30*  loading  step 
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Figure  IS  (b)  The  deformed  mesh  at  the  23^  loading  step 
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Force 


figure  IS  (c)  Hie  displacement  and  the  total  corresponding  nodal  force  curve 


Figure  17  (a)  The  finite  element  mesh  with  the  microcnck  (merited  element) 


(c)  The  deformed  mesh  at  the  29th  loading  «tq> 
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Figure  18  (d)  The  displacement  and  the  total  corresponding  nodal  force  curve 


Figure  19  The  displacement  and  the  total  corresponding  nodal  force  curve 
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Figure  21  Hie  finite  element  mesh  (based  chi  the  Voronoi  diagram)  and  boundary 
conditions 


Figure  22  (c)  The  deformed  mesh  at  the  22nt*  loading  step;  (d)  the  deformed  mesh  at  the  30*  loading  step 


volumetric  Starin 


Figure  24  (a)  The 


Figure  24  (c)  The  deformed  mesh  at  the  25*  loading  step;  (d)  the  deformed  mesh  at  the  30*  loading  step 


Figure  25  The  displacement  and  the  total  corresponding  nodal  force  curves 
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Figure  26  Tbc  grain  size  distribution  using  different  sintering  techniques 
(•)  Normal  sintering,  (b)  Fast  sintering. 


(c)  Long  time  sintering, 
(e)  Two  stage  sintering, 


(d)  High  temperature  sintering, 
(f)  MgO-doped  AI2O3  sintenng. 
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Figure  27  The  void  size  distribution  using  different  sintering  techniques: 

(a)  Normal  sintering,  (b)  Fast  sintering, 

(c)  Long  time  sintering,  (d)  High  temperature  sintering, 

(e)  Two  stage  sintering,  (f)  MgO-doped  AI2O3  sintering. 
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Appendix  I 


The  processing  end  properties  of  alumina 


An  alumina  bar  measuring  1.15  in.  X  1.15  in.  X  1  ft  was  purchased  from  die 
Coon  Technical  Ceramics  Co.  Alumina's  relative  density  is  99.5%. 

The  aluminum  oxide  (or  alumina)  is  prepared  by  hot-pressing  fine  powder  (nominal 

particle  size  OJ  mm),  without  additives,  at  16S0PC  for  3  hours  at  35  MPa  under  a  vacuum. 
The  pressed  material  is  non  porous  (>  99.9%  density),  and  the  final  microstnicture  of 
alumina  is  equiaxed  with  different  grain  sizes  from  5  to  50mm,  which  depends  on  the 
sintering  temperature,  time,  and  pressure.  The  specimens  are  ground  to  the  desired 
thickness  and  die  prospective  viewing  surface  is  polished  with  diamond  paste. 

The  properties  of  die  alumina  bar  are  listed  below  (provided  by  Coors  Technical 
Ceramics  Co.): 


Properties 

AD-99.5  (99.5%  A12Q^ 

Density 

3.88  (g/cm^) 

Surface  finish 

0.9  (mm) 

Grain  size  range 

5  -  50  (mm) 

Water  absorption 

0 

Gas  perm. 

0 

Color 

Ivory 

Flexural  strength 

(MOR) 

20°C 

379  MPa 

Elastic  modulus 

20°C 

372  GPa 

Shear  modulus 

20°C 

152  GPa 

Poisson's  ratio 

20°C 

0.22 

Compressive  strength 

20°C 

2620  MPa 

Hardness 

14.1  GPa 

Tensile  Strength 

25°C 

262  MPa 

Fracture  Toughness 

4.5MPam1/2 
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Appendix  2 


Dynamic  Relaxation 


Nonlinear  structural  analyses  include  materially  nonlinear  problems  such  as 
nonlinear  constitutive  equations  with  small  deformations,  geometrically  nonlinear 
problems  normally  associated  with  buckling  or  a  combination  of  both  types  of 
nonlinearities  (Chen,  and  Schreyer,  1990).  The  nonlinear  analysis  in  the  Unite  dement 
method  can  be  expressed  as 


[K]  (UJ  -  {  F  } 


(1) 


where  [K]  is  the  structure  stiffness  matrix, 

{UJ  is  the  vector  of  nodal  degree  of  freedom, 

(F)  is  the  vector  of  nodal  loads. 

[K]  and  {F}  are  regarded  as  dependent  on  {U }.  The  schemes  for  nonlinear  problems  are 
based  on  step-by-step  load  incrementation  and  an  iteration  procedure  to  correct  the 
linearization. 

The  equation  governing  the  structural  behavior  is  considered  to  be  in  the  form 

<F«J))int  *  (Flext  (2) 

where  {FfU)}^  is  the  vector  of  internal  nodal  force,  which  depends  on  {UJ,  and  {F)ext 
is  the  vector  of  the  external  (applied)  nodal  fence. 

Equation  (2)  is  another  expression  of  Equation  (1),  i.e.,  {F(U)}jnt  =  [K]{UJ,  and 
(F)«tf<F|. 

Since  the  DR  method  is  based  on  structural  dynamic  response,  the  governing 

equation  is  the  appropriate  equation  for  developing  the  DR  method.  For  the  n1*1  time 
increment  the  equation  is  given  by 


[M]  (U)  n+  tq  (U)n  +  [K]  (U)n  -  (F(l”)) 


ext 


i.e.. 


[M]  (ii)n+  (CHU)n  +  {F(Un)}im-{F(tn)}cxt 


(3) 


(4) 


where  [M]  is  the  mass  matrix,  C  is  the  damping  matrix,  t  is  time,  n  indicates  the  n*  time 
increment,  a  superimposed  dot  indicates  a  temporal  derivative,  and  other  terms  are  as 
previously  defined.  To  obtain  the  DR  algorithm,  the  following  central  difference 
expressions  are  used  for  die  temporal  derivatives: 


{U}n’1/2  *(-(U}n'1  +  (U)n)/h 
(U)n  =(-  {U}n*1/2+  {U)n+1/2)/h 


(5) 
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where  h  is  •  fixed  time  increment  The  expression  for  {0} n  is  obtained  by  the  avenge 
value: 

(U)n  «0.5(  (U) n* 1/2 +{0} 0+1/2 )  (6) 

Substituting  Eq.  (5)  and  (6)  into  (4)  gives  the  pair  of  equations  used  to  obtain  next  time 
velocity  and  displacement  as 

(U]n+1  -  (U)n  +  h  (U)n+1/2  (7) 

where  (F)"x[  -  {«!"))„,,  end  |F) 

To  preserve  the  explicit  form  of  Eq.  (7),  [M]  must  be  diagonal  and  to  obtain  the 
form  used  for  DR,  [C]  has  the  form 

[C]  *  c  [M]  (8) 

Substituting  Eq.  (8)  into  (7)  gives 

(C|n+1/2  =  §SJo  (Uln'1/2  +  a>tM]'1((F|"x|  -{F)?t)/(2  +  ch), 

(U)n+1  -lui"  +  b  {U)n+1/2  (9) 

where  [M]'1  is  the  inverse  of  [M].  Since  [M]  is  diagonal,  Eq.(9)  is  algebraic.  That  is,  each 
solution  vector  component  may  be  computed  individually  from 

0  Tm  *  1 u  Tm  *  »  i  •  fJ,,  i )  /  [«a  a*ch)l . 


u?+1  -u?  +h  u|,+1/2 


(10) 


where  the  subscript  i  indicates  the  i^1  vector  component  and  is  the  i*  diagonal  element 
of  [M]. 

In  order  to  start  the  integration,  the  velocity  at  t* 1/2 ,  and  t°  must  be  known.  For 
the  DR  algorithm,  the  starting  conditions  are  of  the  form 


(U)°  #  0;  (U)°«0 
Using  Eq.  (6)  and  the  second  of  Eq.  (1 1)  gives 


(ID 
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(12) 


n)r1/2«-  {ui1/2 


For  die  first  time  increment,  the  first  expression  of  Eq.  (9)  becomes 
(0) 1/2  - h  [Mr1  «F)°„  -  (Fl^, )  /  2 


(13) 


The  central  difference  time  integrator  is  therefore  given  by 


if  n*0; 
if  n  *0; 
for  all  n: 


(U) 


1/2  -hlM]'1  ((F|°x[  -  (F)^  )/2 


(U)"+l/2-§^  |U)n-1/2  +  2h[Mrl([F)“xt  -(Fl^J/C+ch), 
(U)”'fI  -{U)“  +  h  (U)"+1/2  (14) 


To  obtain  the  static  solution  from  the  transient  response  equation,  the  damping 
coefficient,  c,  the  time  increment,  h,  and  the  mass  matrix,  [M],  are  selected  to  obtain  the 
fastest  convergence.  Note  that  only  (F}ext  and  {F}jnt  must  represent  the  physical 

problem,  and  c  and  [M]  need  not  represent  the  physical  structure.  Also,  h  is  dependent  on 
[M]  (Underwood,  1983). 


Formally,  the  DR  algorithm  may  be  written  as: 

(a)  choose  v  ( v  *  ch )  and  [M];  {U)°  given 

(b)  residual  {r)n  =  {F}next  -  {F)nim  , 

(c)  if  {r)n  «  {0}  stop,  otherwise  continue, 

(d)  n  *0;  {UJ 1/2  »  h  [M]'1  (r)°/2  , 
n#0; 


IU}°*0. 


{U}n+1/2-g^  {UJ^  +  ZhtMT1  (r)°/(2  +  v) , 
n+1  *{U}n  +h  (U}n+1/2 


(0  (U) 

(f)  n«n+l;  return  to  (b). 

The  difference  between  the  central  difference  time  integrator  and  the  DR  method  is 
that  v  and  [M]  are  fictitious  values  chosen  so  that  the  static  solution  {r}  «  (0)  is  obtained 
in  a  minimum  number  of  steps.  Also  h  is  a  pseudo* time  increment  which  must  be  chosen  to 
ensure  stability  and  accuracy  of  the  iterations. 

The  internal  nodal  force  {F(U))^t  in  an  element  is: 
where  [KJ c*  fv  |B]^[£][B]dV ,  and  the  superscript  e  defines  the  element 


Matrix  [B]  relates  the  strain  aud  nodal  displacements: 
[el  -  [B]  (UJe 
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where  [e]  is  the  total  strain  vector. 


So.  (F<U))?M.J,[B]T[EHB]dV(U)e 
r 

-  vtB]T[E][BJ(U}edV 

k 

-JvIB]TtE)le]dV 

-pB]T[i)dV 

where  [s]  is  the  total  stress  vector.  In  the  DR  method,  the  stiffness  matrix  [K)  does  not 
need  to  be  formed  and  stored. 


» 


Appendix  3 

Tension -Compression  Loading  Stage  for  FSEM 


Introduction 

For  real-time  fracture  study,  a  closed-loop  control  tension-compression  loading  stage  inside  a 
Fast-Scanning  Electron  Microscope  (FSEM)  is  required.  The  FSEM  is  basically  a  modified  SEM 
with  modified  tools  for  rapid  scan  rates.  Thus  it  is  able  to  observe  specimen  activity  at  greater  than 
1  kHz  framing  rate,  allowing  for  TV  time  observations.  For  the  FSEM  project,  specimens  of 
ceramic  were  tested  in  the  FSEM  chamber.  Image  scanning  was  performed  with  a  digital  data 
acquisition  system,  and  the  resulting  micrographs  were  later  examined  in  an  observable  time  frame. 
The  major  drawbacks  of  FSEM  microscopy  include  sacrificed  depth-of-field  and  submicron 
resolution,  resulting  in  a  snowier  picture  than  conventional  SEM  micrographs. 

Objectives 

The  following  list  outlines  the  objectives  of  this  project. 

1.  Assess  use  of  pre-existing  equipment  and  ideas 

2.  Achieve  open-loop  control  of  die  loading  stage,  complete  with  sensor  data  acquisition 

3.  Advance  the  system  to  a  closed-loop  status  by  using  a  selected  motion  control  sequence 

4.  Install  the  completed  in-si tu  device  for  use  inside  the  SEM  chamber 

5.  Perform  tests  with  the  alumina  ceramic  specimen  to  observe  effects  of  the  work 

General  Principle  and  System  Elements 
General  Principle 

The  whole  experimental  system  is  controlled  by  die  Galil  Motion  Controller  DMC-1010  which 
plugs  into  the  PC  bus.  It  sends  designated  commands  to  the  step-motor  set,  which  drives  the 
worm-gear  assembly  and  screw  shaft  to  move  the  specimen  platform.  While  die  specimen  is  being 
dragged,  a  clip-on  gage,  which  is  held  between  knife  edges  attached  to  the  specimen,  moves 
simultaneously  to  sense  the  crack  opening  displacements.  A  load  cell  is  also  engaged  to  test  the 
changing  tension  on  the  specimen.  The  data  from  the  click-on  gage  and  load  cell  are  amplified  and 
filtered  and  then  collected  by  tire  motion  controller.  In  a  closed-loop  control  test,  the  sample  data 
are  analyzed  by  the  controller  and  an  error  signal  is  sent  out  to  adjust  the  system's  movement  The 
Visual  Basic  system  is  introduced  to  establish  the  testing  data  file  and  to  paternally  control  the 
testing  system.  The  working  system  is  shown  in  Figure  3-1. 


System  Elements 

The  Galil  Motion  Controller  DMC-1010  is  designed  for  maximum  system  flexibility,  and  is 
available  for  one,  two,  three,  or  four  axes.  It  can  be  interfaced  to  a  variety  of  motors  and  drives 
including  full-step,  micro-stepping  and  servo-motor.  Its  main  processing  unit  is  a  specialized  32- 
bit  Motorola  68331  Serbs  Microcomputer  with  64  K  RAM,  64  K  EPROM  and  256  bytes  EPROM. 
It  provides  interface  circuitry  for  8  outputs,  8  inputs  and  7  analog  inputs.  The  analog  input  is  a  12- 
bit  ADC  with  a  range  of  ±  10  V,  which  is  0.488  mV/bit. 

The  step-motor  set,  which  is  an  American  Precision  Industries  model  P261-M232,  includes  a 
drive  controller  and  a  step  motor.  It  combines  an  integral  power  section  and  the  bipolar  chopper 
drive  technique.  This  technique  involves  over-driving  the  windings  with  high  voltages  to  decrease 
the  current  by  high  frequency  current  chopping.  The  number  of  steps  per  revolution  of  the  motor 
can  be  set  with  dip  switches  on  the  controller.  The  current  setting  is  50800  step/revolution. 

Controlling  the  motor  requires  three  basic  inputs  to  the  controller.  The  first  input  is  the  step 
input.  The  motor  will  step  on  the  trailing  edge  of  each  incoming  pulse  up  to  a  rate  of  1  MHz.  The 
pulse  width  is  set  in  DMC-1010  software  to  1920  NS.  The  second  input  directs  motor  rotation.  In 
this  test,  voltage  is  taken  low  corresponding  to  the  counterclockwise  (CCW)  rotation.  The  third 
input  to  the  controller  is  an  OPTP  input  which  requires  +  5  VDC  to  operate  the  optical  isolation 
feature.  These  three  inputs  are  controlled  and  provided  by  DMC-1010. 

The  specimen  test  set  includes  the  worm-gear  assembly,  the  specimen  mounting  stage,  the 
specimen  platforms,  the  clip-on  displacement  gages  and  the  load  cell.  The  platforms  are  moved 
inward  or  outward  by  the  worm-gear  and  the  treaded  screw.  The  screw  and  the  platforms  are 
mounted  on  the  stage.  The  specimen  is  30  mm  square  and  1  to  3  mm  thick  with  two  3-mm 
diameter  holes  drilled  about  14  mm  apart  to  fit  over  the  loading  pins  which  are  screwed  into  the 
platforms  (see  Figures  3-2  and  3-3). 

The  Clip-On  Displacement  (COD)  Gage  is  a  product  of  MTS  Systems  Corporation  (MTS 
model  632.03B-30).  It  is  designed  to  sense  crack  opening  displacements  and  is  typically  used  to 
monitor  crack  growth  in  a  specimen  during  fracture  mechanics  testing.  It  uses  precision  resistance- 
type  strain  gages  bonded  to  a  metallic  element  to  form  a  whetstone  bridge  circuit  The  calibration 
of  COD  is  0.007131  in  ./mV,  or  0.18  mm/mV,  on  a  4  mm  scale.  The  excitation  voltage  is  5.99  V. 

The  Load  Cell  (LC)  is  a  product  of  Sensotec,  Inc.,  part  no.  1 1/2335-07  and  serial  no.  353700. 
The  LC's  full-scale  capacity  is  1000  lbs  under  tension.  The  calibration  factor  is  2.4602  mV/V. 
Because  the  excitation  is  5.0  V,  the  calibration  is  2.4602  mV/V  x  5.0  V  *  12.3  V  for  1000  lbs,  or 
81.3  lbs/mV. 

Hie  amplifiers  are  transducer  conditioning  amplifiers  manufactured  by  die  Ectron  Corporation. 
The  two  Model  563  FN  amplifiers  are  inducted  in  the  Model  E5 13-20  Enclosure.  They  are  wide 
band,  true  differential  DC  instrumentation  amplifiers  with  built-in  transducer  signal  conditioning 


functions.  The  output  of  the  basic  amplifier  is  ±  10  V  at  10  mA.  The  maximum  gain  is  1  K.  The 
amplifiers  used  in  this  test  were  for  COD  and  LC  signal  amplifying  and  for  the  supply  of  excitation 
power. 

The  filter  is  a  Krohn-Hite  Model  3988  Butterworth/Bessel  dual  channel  filter.  It  provides  a 
tunable  frequency  range  form  0.03  Hz  to  1  MHz  in  low-pass  mode  and  0.03  Hz  to  300  kHz  in 
high-pass  mode.  Each  channel  of  the  3988  is  an  8-pole,  wide  range,  low-pass/high-pass  filter  or 
an  amplifier  providing  gains  to  70  dB  in  0. 1  dB  steps.  The  3988  will  accept  input  signals  of  ±  10 
V  peak  at  0  dB  gain  and  has  selectable  ac  or  dc  coupling.  The  two  channels  of  the  filter  are  used 
for  COD  and  LC  signals  in  this  test 

System  Noise  and  Noise  Solving 

Noise  is  the  main  problem  in  this  micro-displacement  test  The  noise  in  the  system  is  more 
distinguishable  than  the  useful  signals.  Three  kinds  of  noise  were  found  to  exist  using  an 
oscilloscope  and  FFT  equipment  to  analyze  the  system  noise. 

One  type  of  noise  was  caused  by  the  discrimination  of  the  Galil  Controller.  The  lowest  limit 
for  the  ADC  of  DMC-1010  is  0.488  mV;  therefore,  the  sampling  signals  should  be  much  higher. 
This  problem  can  be  effectively  eliminated  by  amplifying  the  useful  signals  before  they  are  sampled 
by  the  ADC.  Considering  the  highest  limit  (±  10  V)  of  the  ADC,  the  gain  for  the  displacement 
signals  was  set  to  500  and  the  gain  for  load  signals  to  1000.  Thereafter,  the  useful  signals  were  in 
the  range  of  voltage  and  the  ADC  noise  was  suppressed. 

The  second  kind  of  noise  was  caused  by  electricity.  The  sources  of  this  noise-causing 
electricity  included  the  computer,  the  step-motor,  the  motor  controller,  and  the  power  supply 
source.  By  using  FFT  equipment,  two  powerful  noise  signals  were  found.  One  signal  was 
around  60  Hz  and  the  other  was  around  1.5  kHz.  They  were  tested  when  the  step  motor  was 
turned  off. 

When  the  step  motor  was  on,  the  third  kind  of  noise  appeared.  The  strongest  AC  noise  signal 
was  12.5  Hz  with  a  magnitude  of  0.076  mV.  Another  AC  noise  signal  was  around  1  Hz,  or  0.016 
mV.  These  noise  signals  were  tested  without  amplifying.  They  were  probably  mechanical  noise. 

An  attempt  was  made  to  eliminate  the  noise  by  shielding  the  step-motor  controller  bar  and 
covering  the  motor  set  and  line  with  tin  foil.  The  best  results  were  obtained  after  using  a  low-pass 
filter.  The  filter  was  set  for  a  low-pass  filter  of  0.5  Hz  and  for  DC  signals.  Although  a  lower 
frequency  noise  signal  exists  with  a  maximum  of  0.02  mV,  0.5  Hz  is  still  an  ideal  cutoff  frequency 
because  any  lower  filter  cutoff  frequency  causes  the  system  response  to  delay  severely  in  open- 
loop  and  closed-loop  control  tests. 


Open-Loop  Control 

In  the  open-loop  control  test,  specimens  were  made  using  both  glass-plastic  materials  and 
aluminum  ceramic  material.  Because  the  glass-plastic  specimens  are  easier  to  make,  they  were 
used  in  the  first  period  of  test  In  the  test  the  step  motor  was  set  to  rotate  at  3500  counts  per 
second.  At  that  rotation  rate,  the  specimen  was  dragged  by  a  speed  of  0.002  mm/sec,  which  is 
0.01 12  mV/s  to  the  COD  correspondingly.  The  sampling  period  was  decided  by  the  characteristics 
of  glass-plastic  and  the  limitation  of  the  DMC-1010. 

The  glass-plastic  had  a  relatively  large  strain.  The  size  of  the  specimen  was  30  x  30  x  2.5  mm. 
It  was  made  in  the  same  way  as  the  ceramic  specimen  was  made.  It  would  not  break  until  the  total 
strain  was  above  0.25  mm.  The  DMC-1010  supplied  an  array  for  a  total  of  1600  elements. 
Therefore,  if  the  COD  signals  and  LC  signals  were  sampled  simultaneously,  a  maximum  of  800 
data  samples  could  be  collected  for  each.  To  record  the  data  of  all  the  periods  of  a  test,  1  second 
was  chosen  as  the  sampling  period.  The  load-displacement  graph  is  shown  in  Figure  3-4.  The 
alumina  ceramic  material  was  more  brittle  than  the  glass-plastic  material.  When  the  specimens 
measured  30  x  30  mm  with  a  thickness  of  0.9-1.96,  they  sustained  a  strain  amount  of  about 
0.04-0.06  mm  in  a  range  of  0.2208-0.3324  mV  of  COD. 

Several  step  motor  speeds  were  chosen  for  the  tests.  One  speed  was  designated  to  be  500 
counts/second.  Thus  the  specimens  were  loaded  at  a  constant  load-point  displacement  of  0.43  pm 
per  second,  corresponding  to  0.0024  mV/s  of  COD.  Another  speed  was  set  at  1000  counts/sec,  or 
0.0048  mV/s.  The  load-displacement  graphs  are  shown  in  Figures  3-5  and  3-6.  From  the  load- 
displacement  curve,  it  can  be  observed  that  the  higher  the  step  motor  speed,  the  more  linear  the 
curve.  This  feature  occurred  because  the  stiffness  of  the  test  system  was  not  hard  enough, 
especially  the  two  pins  which  held  die  specimens.  Since  the  speed  of  the  test  is  still  relatively  high, 
it  is  difficult  to  obtain  an  ideal  characteristic  curve  of  load-displacement  for  alumina  ceramic 
material. 

Closed-Loop  Control 

The  principle  of  closed-loop  control  is  shown  in  Figure  3-1.  The  output  signals  were  obtained 
from  the  COD  and  LC,  and  then  the  error  signals  were  fed  back  into  the  system.  Because  the 
rotating  speed  of  the  step-motor  was  constant  in  the  test,  the  distance  between  two  sample  points 
was  chosen  for  reference.  On  this  basis,  the  average  velocity  between  every  two  sample  points 
was  assumed  to  be  the  same.  If  there  was  any  difference,  the  error  was  fed  back  into  the  step- 
motor  to  increase  or  decrease  the  rotation  speed  of  the  motor.  For  the  glass-plastic  specimens, 
when  the  motor  rotated  at  a  speed  of  3500  c/s,  the  change  of  COD  was  0.01 12  mV/s  as  mentioned 
before.  Because  the  noise  was  about  0.02  mV,  the  sample  period  was  set  to  be  2.5  seconds.  Thus 
the  displacement  changes  were  0.028  mV  for  COD  in  each  sampling  interval.  Then  the  range  to  be 
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controlled  was  from  0.028  -  0.01  mV  to  0.028  +  0.01  mV.  Figure  3-7  shows  the  load- 
displacement  curve  used  in  this  control  range.  Compared  with  the  results  from  the  open-loop  test, 
the  maximum  tension  and  maximum  strain  displacement  are  almost  the  same;  however,  the  curve 
from  the  closed-loop  control  system  is  much  more  linear. 

Difficulties  arise  when  using  the  existing  test  equipment  to  perform  closed-loop  control  tests  on 
ceramic  samples.  As  analyzed  previously,  the  maximum  strain  of  ceramic  is  about  0.04-0.06  mm, 
or  0.2208-0.3324  of  COD,  but  the  noise  signal  is  0.02  mV.  Therefore,  the  noise  signal 
overwhelms  the  useful  signal  in  micro-displacement  movement 

Conclusion 

The  system  is  suitable  to  perform  open-loop  tests  for  the  micro-displacement  test  The  noise  in 
the  system  is  still  a  significant  problem.  It  may  be  caused  by  the  mechanical  equipment.  The  test 
system  has  difficulty  performing  more  precise  testing,  especially  closed-loop  control  tests. 
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Step  Motor  Set 


Figure  £-1.  System  Block  Diagram 
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S  -  Specimen 

G  -  Worm  Gear 

P  “  Loading  pins 

C  -  Load  cell 

LP  -  Loading  platforms 

B  -  Base  of  loading  stage 


Figure  3-2.  Specimen  Test  Set 
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figure  3-3.  Thread  Screw  Bar 
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Figure  3-6.  Open-Loop  Control 

load- Displacement  for  Alumina  Ceramic  Specimen 
Motor  rotation  speed  1 000  counts/second 
Samp&ng  period  1  seconds,  LPF  0.S  Hz 
Sample  thickness  1.98  mm 
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Figure  3-7.  Close- Loop  Control 

Load-Despla cement  Graph  for  Glass-Plastic  Specimen 


APPBMDIX  4 


Product  Specification 

Timothy  J.  Ross  and  Ming  L.  Wang,  University  of  New  Mexico 


"FSEM:  FAST  SCANNING  ELECTRON  MICROSCOPY" 

U.S.  Pat.  No.  5,254,857  Oct.  1993 

Why  High-Speed  SEM? 

Transmission,  reflection  and  scanning  electron  microscopes  offer  spatial  resolution 
superior  to  that  of  light  microscopes.  Transmission  and  reflection  electron  microscopes  also 
provide  temporal  resolution  of  -20  ns/frame,  rivaling  the  temporal  resolution  of  light  microscopes. 
SEM,  however,  is  the  only  type  of  electron  microscope  which  offers  both  high  spatial  resolution 
and  large  depth  of  focus  surface  imaging,  even  for  rough  surfaces,  of  particular  value  for  surface 
studies  of  material  responses  to  dynamic  loads  at  high  magnifications,  for  example.  Based  on 
experiments  with  the  high-imaging  speed  (FSEM:  Fast  SEM)  now  under  development  at  the 
University  of  New  Mexico  (UNM)  under  a  collaboration  with  industry,  we  export  that  state-of-the- 
art  components  will  soon  permit  SEM  dynamic  microscopy  at  0.1  pm  spatial  resolution  and  2 
p s/frame  at  64  pixels  X  64  pixels  per  frame,  with  acceptable  signal-to-noise  ratios  in  final  images. 

Problems  of  High-Imaeing  Speed  SEM 

To  obtain  high-imaging  speed  operation  in  the  SEM,  four  task  areas  must  be  addressed: 
illumination,  detection,  deflection  and  recording.  The  specimen  must  be  illuminated  with  an 
electron  beam  of  sufficient  intensity  to  provide  a  detectable  signal  of  sufficient  quality  to  provide  a 
good  image  of  the  specimen.  And,  such  an  image  must  be  obtained  at  framing  speeds  of  interest 
with  scan  beam  deflection  providing  suitable  speed  and  proper  deflection  waveform,  without 
degrading  scan  beam  quality.  The  fourth  task  is  to  record  the  imaging  signal  synchronously  with 
the  scan  for  later  playback  in  a  "movie"  format.  To  address  these  four  problems  and  to  illustrate 
the  FSEM  technology,  we  modified  an  SX-40A  SEM  manufactured  by  International  Scientific 
Instruments,  Inc.  The  unmodified  SEM  had  a  tungsten  hairpin  cathode  electron  gun,  two  stages  of 
condenser  optics,  two-stage  magnetic  deflection  and  an  objective  lens. 

Illumination  and  Detection 

Illumination  and  detection  are  intimately  related.  If  the  detector  detects  all  imaging  radiation 
produced  by  the  scan  beam  without  introducing  noise  or  frequency  errors,  funher  image  quality 
improvements  can  only  be  had  by  increasing  illumination.  Though  our  detector  is  not  this  perfect, 
we  believe  its  deficiencies  have  far  less  impact  On  image  quality  than  those  of  the  illumination 
system  we  have  at  present.  In  any  case,  illumination  ultimately  determines  obtainable  image 
quality.  In  the  case  of  secondary  electron  imaging,  for  example,  and  typical  scan  beam  (primary 
electron)  energies,  one  secondary  electron  is  produced  for  every  1-10  primary  electrons.  Clearly, 
improving  primary  beam  current  can  dramatically  improve  image  quality. 
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Deflection 


The  other  major  problem  area  we  have  encountered  in  obtaining  high-imaging  speed  SEM 
operation  is  in  deflection  of  the  scan  beam  at  high  rates  with  minimal  distortion  of  the  desired 
deflection  waveform  (non-linearity)  or  degradation  erf'  scan  beam  quality  (deflection  aberrations). 
Several  improvements  to  scan  generator/scan  coil  driver  electronics  were  performed  Problems  due 
to  the  quiescent  state  of  the  scan  generator  producing  a  constant  DC  input  to  the  scan  coil  drivers 
between  scans,  resulting  in  image  shift  and  unnecessary  power  dissipation  in  the  scan  coil  drivers 
and  scan  coil  terminadon  resistors  have  been  resolved.  Lower  magnification  was  required  to 
properly  center  specimens  in  the  microscope  field  of  view  for  dynamic  imaging,  but  the  lower 
magnification  requires  higher  deflection  currents  causing  significant  heating  in  the  scan  coil 
termination  resistors.  These  problems  were  also  resolve!  Another  modification  which  makes  the 
microscope  easier  to  use  is  a  switch  box  permitting  simple  switching  between  analog  data 
acquisition  and  digital  data  acquisition.  The  specimen  is  first  positioned  in  the  microscope  field  of 
view  in  analog  mode,  then  switched  to  digital  mode  to  take  a  digital  movie. 

Rscoiding 

The  imaging  signal  output  is  input  to  a  digital  data  acquisition  system  ,  which  consists  of  a 
fast  digitizer,  a  digitizer  controller  (via  GPIB  bus),  fast  memory  modules  and  a  chassis  containing 
appropriate  power  supplies.  The  digitizer  we  are  using  at  present  is  the  LeCroy  8828D,  capable  of 
a  200  MHz  digitizing  rate,  8-bit  measurement  and  1  Mbyte  of  fast  memory  storage  (expandable  to 
2  Mbyte).  At  256  pixels  X  256  pixels  per  frame  this  allows  for  storage  of  16  frames  of  data.  At 
128  X  128  pixels  per  frame,  64  frames  of  data  may  be  stored.  The  digitizer  also  provides  a  clock 
pulse  used  to  generate  scan  deflection  signals  so  that  data-taking  and  scanning  are  synchronized. 

Project  Status 

In  the  last  year  we  have  obtained  images  at  a  25  MHz  pixel  rate,  which  translates,  for  256 
pixels  X  128  pixels  per  frame,  to  381  frames  per  second.  This  is  -10  times  TV-rate  (nominally 
30-50  Hz),  the  maximum  framing  rate  available  in  any  commercial  SEM.  Further,  to  our 
knowledge,  25  MHz  pixel  rate  is  -2-3  times  faster  than  anyone  else  has  obtained  for  secondary 
electron  or  any  other  sort  of  SEM  imaging.  We  have  developed  a  design  which  will  allow  us  to 
obtain  images  at  the  200  MHz  pixel  rate  specification  of  our  digital  data  acquisition  system  and  at 
higher  rates  using  a  faster  digital  data  acquisition  system.  A  200  MHz  pixel  rate  translates  to  a 
12.2  kHz  framing  rate  at  128  pixels  X  128  pixels  per  frame.  A  2  GHz  pixel  rate  (state-of-the-art  in 
8-bit,  long  record  length  digitizers)  results  in  a  488-kHz  framing  rate  at  64  pixels  X  64  pixels.  Our 
images  have  acceptable  signal-to-noise  ratios,  i.e.  S/N  >  6. 

We  have  obtained  images  at  a  381  Hz  flaming  rate,  a  beam  voltage  of  10  kV  and 
magnification  of  >  1.000X.  We  believe,  however,  that  we  are  at  the  limit  of  the  image  quality  that 
can  be  obtained  with  the  electron  source  we  have  at  present,  and  at  the  stated  operating  parameters. 
Some  brightness  improvement  is  possible  by  choosing  different  operating  parameters.  However, 
by  using  high  brightness  cathodes,  which  may  be  2-3  orders  of  magnitude  brighter  than  the  source 
we  have  at  present,  we  believe  that  image  quality  may  be  greatly  improved.  We  have  also  obtained 
time-integraied  photographs  taken  at  various  horizontal  line  times  u>  examine  high-speed  linearity 
and  aberrations  of  the  deflection  coils,  as  well  as  the  time  resolution  of  the  high-speed  secondary 

electron  detector.  If  we  assume  128  pixels  per  horizontal  line,  a  horizontal  line  time  of  7.5  ps 
represents  a  time  resolution  of  60  ns.  Also  assuming  128  vertical  pixels  per  frame,  the 
corresponding  flaming  rate  would  then  be  - 1  kHz.  Table  1  shows  the  relationship  of  effective 
framing  rates  for  various  horizontal  line  times. 


With  the  FSEM  capabilities  discussed  above,  a  variety  of  dynamic  microscopic  studies 
become  possible.  Chief  among  these  is  microstructural  response  of  structural  materials  to  transient 
loads,  such  as  mechanical,  thermal  and  magnetic  loads.  We  have  developed  a  magnetically- 
induced  stress  wave  (MISW)  device  for  producing  mechanical  loads  in  specimens.  This  device 
has  fractured  DSP  concrete  cement  in  dynamic  tension  loading.  The  pressure  pulse-width  of  this 

device  is  ~1  jis,  with  pressure  pulse  peaks  of  £  10,000  psi.  This  magnetically-induced  stress 
wave  concept  using  -1  mm  radius  rod  specimens  allows  for  various  combinations  of  dynamic 
loading-single-  or  double-ended,  compression  or  tension-- and  does  permit  fracture  to  be  localized 
within  the  SEM  field  of  view  and  to  time  fracture  with  scanning.  However,  these  advantages  can 
only  be  realized  for  specimen  materials  which  are  homogeneous  (on  the  scale  of  the  pressure  pulse 

length)  and  (i)  can  transmit  a  ~1  ns  wide  pressure  pulse  without  dissipation,  (ii)  are  weak  enough 

to  break  at  pressure  peaks  in  die  3,000-30,000  psi  range  (the  present  limits  of-our  device - 

extrapolated  to  the  35  kV  maximum  capacitor  charge  voltage)  and  (iii)  exhibit  crack  nucleation  and 

growth  and  macroscopic  fracture  within  -1  (is.  The  FSEM  imaging  is  synchronized  with  the 
initiation  of  fracture  in  the  specimen. 
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During  the  research  effort,  the  concept  of  using  a  high-speed  scanning  electron  microscope 
(SEM)  to  observe  real-time  microstmctural  response  of  dynamically  loaded  structural  materials  was 
verified  experimentally  at  a  maximum  framing  rate  of  381  Hz  (256  horizontal  pixels  X  128  vertical 
pixels),  about  an  order  of  magnitude  higher  than  previously  possible  with  conventional  SEM's. 
This  experimental  accomplishment  proved  the  soundness  of  several  key  concepts: 

*  That  a  tungsten  hairpin  cathode  is  bright  enough  to  obtain  useful  digital  images  at 
the  framing  rate  listed  above, 

*  That  a  secondary  electron  detector  can  be  built  and  operated  at  high  enough  count 
rates  to  obtain  such  images, 

*  That  the  scan  coil  assembly  standard  on  an  IS  I SX-40A  SEM  can  be  replaced  to 
allow  imaging  at  such  rates  with  spatial  resolution  approaching  100  nm, 

*  That  signal  acquisition  and  scan  generation  can  be  synchronized  to  obtain  a 
succession  of  well-defined  frames  in  a  “movie"  format  at  pixel  rates  far  in  excess  of 
conventional  TV -rate  SEM  video  band  widths  (Table  2  lists  the  events  for  which  we 
have  made  movies),  and 

*  That  a  magnetically-induced  stress  wave  device  can  be  used  to  obtain  dynamic 
fracture  within  the  SEM  chamber  and  field  of  view,  with  scanning  timed  to  coincide 
with  fracture. 

Based  on  the  progress  made  during  this  research  effort  and  the  promise  that  the  FSEM  device  has 
for  the  new  area  of  dynamic  microscopy  a  patent  has  been  filed  on  the  technology. 


•  Glass  beads  excited  by  piezoelectric  transducer 

•  Thomel  fiber  excited  by  piezoelectric  transducer 

•  Fuse  rupture  under  excess  current  pulse 

•  Rosin  fracture  under  shock  loading 

•  Concrete  fracture  under  tensile  load 

•  Watch  gear  moving  periodicity 


Table  3  shows  the  technology  areas  that  we  have  reviewed  over  the  last  2  years  that  could 
benefit  from  die  utility  of  our  FSEM  device.  Certain  modifications,  such  as  thermal  stages  or  a 
gaseous  environment,  to  the  stages  would  be  necessary,  but  the  FSEM  technology  could  be  readily 
adapted  to  these  situations. 


>  monolayer  studies  (Langmuir-Blodgett)*-may  require  ESEM 
1  recrystallization  of  silicon  films  and  semiconductors 
1  real-time  dynamic  fracture 
1  In  situ  SEM-  or  microprobe-observed  fracture 
■  TEM  video  microscopy  (TV-rate) 

1  solid  phase  transitions  (see  high  Tc  superconductor  listings  below) 
atomic  motion  within  crystals 

Pulsed  electron  microscopy-non-scanned  (—25  ns/frame) 

TEM  studies  of  laser-pulse-induced  processes 
reflection  electron  microscope  surface  studies 
High  Tc  superconductor  studies 
phase  transitions  during  fabrication-SEM 
flux  creep,  magnetic  contrast  imaging-SEM 
flux  creep,  decoration  imaging-SEM 
Domain  growth/transition  of  polymers 
Integrated  circuit  inspection  and  metrology 

E-beam  testing  of  integrated  circuits  (strobed  voltage  contrast  in  SEM) 
E-beam  specimen  damage  reduction 
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ABSTRACT 


Constraints  arise  naturally  in  the  context  of  rigid  inclusions, 
incompressibility,  inextensible  fibers  and  deformed  finite  elements.  Typical 
methods  for  handling  constraints  in  the  governing  matrix  equation  include 
Lagrange  multipliers,  penalty  weights,  and  elimination  of  variables.  Each  has  a 
particular  undesirable  feature.  Proposed  here  is  a  procedure  in  which  each 
constraint  is  handled  directly  and  sequentially  through  a  modification  of  the  rows 
and  columns  of  the  governing  matrix  and  force  vector.  Positive  definiteness, 
symmetry  and  the  dimensions  of  the  matrix  remain  unchanged.  Elementary 
examples  involving  both  the  static  and  dynamic  response  of  a  bar  are  given  to 
illustrate  the  procedure. 


Key  Words:  constraints,  rigid  inclusions,  incompressibility,  dynamic  response 
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INTRODUCTION 


There  are  numerous  engineering  situations  in  which  constraints  are  a 
natural  part  of  the  system.  Examples  are  rigid  inclusions,  incompressibility, 
inextensiole  fibers,  floors  in  plane  frames,  boundary  conditions,  and  mesh 
distortion.  The  usual  methods  for  handling  constraints  require  either  a 
modification  to  the  dimension  of  the  algebraic  problem,  or  good  numerical 
intuition.  On  the  other  hand,  the  implementation  of  a  one-point  constraint  is 
often  done  directly  by  merely  adjusting  rows  and  columns  so  that  the  original 
dimension  and  symmetry  of  the  matnx  remain  unchanged.  A  corresponding 
procedure  for  a  multipoint  constraint  of  arbitrary  dimension  is  not  generally 
available,  at  least  in  most  textbooks  on  the  finite  element  method.  Here,  such  an 
algorithm  is  given  with  simple  examples  for  static  and  dynamic  problems. 

If  multipoint  constraints  are  given  explicitly,  conventional  methods  for 
handling  the  constraints  are  the  Lagrange  multiplier  technique,  the  penalty 
method,  and  partitioning  of  variables  with  the  subsequent  need  for  adjustments 
to  the  governing  matrix  (Barlow,  1982;  Cook  et  al.,  1989).  The  Lagrange- 
multiplier  method  requires  the  development  of  an  augmented  matrix  which  means 
the  problem  must  be  redimensioned.  The  penalty  method  involves  an  adroit 
choice  of  a  parameter  to  provide  a  balance  between  ill  conditioning  and  loss  of 
accuracy.  Partitioning  of  variables  requires  a  considerable  amount  of 
manipulation  and  the  inverse  of  a  submatrix  which  can  be  of  a  significant  size  if 
there  are  a  large  number  of  constraints. 

Sometimes  governing  systems  become  ill-conditioned  as  a  parameter 
describing  a  material  property,  the  dimension  of  a  body,  or  element  distortion 
approaches  a  critical  value.  Various  techniques  involving  under  integration  with 
hourglass  control,  other  modifications  to  the  element  stiffness  matrix,  and 
remeshing  are  often  used  to  prevent  the  ill  conditioning.  For  many  cases,  these 
techniques  can  be  interpreted  as  approximate  methods  for  applying  a  constraint 
which  should  be  applied  based  on  fundamental  mathematical  arguments. 

Webb  (1990)  uses  a  projection  operator  to  elegantly  derive  a  rather 
general  algorithm  for  handling  constraints  without  modifying  the  structure  or  the 
symmetry  of  the  governing  matrix.  However,  nothing  is  said  about  how  the 
projection  is  generated.  Hueck  and  Schreyer  (1993)  provide  an  explicit  method 
for  constructing  the  projection  matrices.  The  result  is  a  general  method  for 
exactly  handling  multipoint  constraints  without  altering  the  number  of  algebraic 
equations.  As  an  example,  the  incompressiblity  constraint  is  invoked  for  the 
quadrilateral  element  with  plane-strain  elasticity.  With  this  approach,  the  problem 
normally  associated  with  the  usual  method  of  fetting  Poisson’s  ratio  approach  0.5 
is  avoided.  Also,  severe  mesh  distortion  is  not  a  significant  factor.  Although  the 
application  to  a  beam  problem  provides  good  results,  the  projection  method 
suffers  from  the  fact  that  a  matrix  as  large  as  the  original  one  and  the  inverse  of  a 
smaller  matrix  have  to  be  obtained.  Therefore,  the  potential  usefulness  of  the 
approach  for  large-scale  problems  is  limited. 

To  maintain  the  advantage  of  the  direct  application  of  multipoint 
constraints,  an  algorithm  is  proposed  that  is  closely  related  to  one  given  by  Abel 
and  Shephard  (1979),  a  contribution  which  has  not  received  the  attention  it 
deserves.  The  attractive  features  of  the  method  presented  here  include  the  fact 
that  the  symmetry  and  the  dimension  of  the  governing  matrix  is  maintained,  there 
is  no  reordering  of  variables,  and  with  some  care,  any  number  of  constraints  can 
be  applied  sequentially  with  the  one  algorithm.  A  complete  description  of  the 
procedure  is  given  together  with  static  and  dynamic  examples. 


BASIC  EQUATIONS 


Suppose  the  objective  is  to  solve  for  x  the  algebraic  equation 

Ax  »  b  (1) 

in  which  the  symmetric  matrix  A  is  n  x  n  and  b  is  given.  The  system  is  restricted  by 
the  constraint 

cTx  «  d  (2) 

An  alternative  way  of  viewing  the  problem  is  to  define  the  residual 

r  «  Ax-b  (3) 

Consider  a  weighting  vector  w.  Frequently  the  governing  equation  (1)  is  obtained 
from  the  requirement  that 

wTr  *  0  (4) 

for  arbitraiy  w.  Whenever  a  component  of  x  is  prescribed,  the  corresponding 
component  of  w  must  be  zero.  Similarly,  if  the  constraint  (2)  exists,  then  w  is  not 
arbitrary  but  must  satisfy  the  corresponding  constraint 

cTw  *  0  (5) 

The  formulation  given  by  (4)  and  (5)  is  the  one  that  provides  insight  into  how  the 
matrix  A  should  be  modified  so  as  to  maintain  symmetry.  If  A  is  not  initially 
symmetric,  then  the  modified  matrix  is  also  not  symmetric. 

For  convenience,  let  {A;}  and  <A,>  denote  the  i’th  column  and  row  of  A, 
respectively.  Then  an  alternative  form  of  (4)  is 


£wiri 

i=l 


Ti  =  <Aj>{x}  -bj 


l  ™  l,...,n 


handle  the  one-point  constraint  normally 
Then  the  general  algorithm  is  given. 


THREE-POINT  CONSTRAINT 


For  ease  of  presentation,  the  constraint  is  applied  to  the  first  three 
components  of  x  but  he  algorithm  can  be  applied  equally  well  to  any  other 
combination  of  components. 

Suppose  the  constraint  is 

CiXi  +  C2X2  +C3X3  ■  d  (7) 


The  constraint  on  w  is 
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ClWj  +  C2W2  +  C3W3  *  0 


(8) 


which  implies  that  wj,  W2  and  W3  are  not  arbitrary.  Set  W3  ■  -C1W1/C3  -  C2W2/Q. 
The  result  is  that  each  of  the  first  two  equations  is  really  a  linear  combination  of 
the  original  equation  and  the  third  equation.  Since  W3  has  been  eliminated,  the 
constraint  equation  is  used  for  the  third  equation.  At  this  stage,  the  governing  set 
of  algebraic  equations  is: 


A*x  -  b' 


(9) 


in  which 
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Now  suppose  the  constraint  is  used  to  replace  X3  in  every  equation  but  the  third. 
Then  the  governing  matix  and  force  vector  become: 
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in  which  a*ij  and  b*j  denote  the  components  of  A*  and  b*,  respectively.  The  only 
components  not  symmetric  involve  tne  indices  (1,  2),  (13)  and  (2,  3).  Suppose 
the  third  equation  is  multiplied  by  C1/C3  and  the  result  added  to  the  first  equation. 
Then  the  tnird  component  in  the  first  row  of  A**  is  ci.  Similarly,  multiply  the 
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third  equation  by  C2/C3  and  add  the  result  to  the  second  equation  to  obtain  C2  as 
the  third  component  in  the  second  row  of  A**.  Now  the  governing  matrices  are 
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It  can  be  shown  by  direct  substitution  that  the  (1,  2)  and  (2,  1)  components  of  the 
modified  matrix  are  equal  so  the  resulting  system  is  symmetric.  The  solution  to 


*•*, 


(16) 


automatically  satisfies  the  three-point  constraint  equation. 


GENERAL  ALGORITHM  FOR  MULTIPOINT  CONSTRAINTS 

• 

The  results  for  the  special  case  given  above  suggests  the  following 
approach  to  handle  a  general  multipoint  constraint.  The  lormat  is  chosen  to 
suggest  the  form  that  a  subroutine  might  take  for  invoking  a  given  constraint.  Let 
Sc  be  the  set  of  indices  associated  with  the  constraint.  Then  tne  constraint  is  given 
by 

£qxj  *  d  (17) 

i«Sc 

in  which  the  summation  is  taken  over  all  indices  in  the  constrained  set.  Each  q  is 
nonzero.  Choose  arbitrarily  a  particular  element  of  S0  i  *  I  say,  with  xj  labelled 

•  the  constraint  variable  {or  master  variable).  One  might  wish  to  choose  I  so  that  q 
is  the  one  with  the  maximum  absolute  value  of  all  q’s  in  the  set  Let  S’ c  be  the  set 
Sc  with  element  I  excluded,  i.e.,  S*c  is  the  set  of  indices  associated  with  the  slave 
variables. 

For  ease  of  interpretation  in  the  three-point  constraint  case,  different 
symbols  were  used  at  eacn  step  for  the  modified  matrices  and  vectors  (e.g..  A*,  b*, 

•  etc.).  However,  there  is  no  need  to  actually  defne  new  matrices  and  vectors  so  in 
the  following  algorithm  only  one  matrix  and  vector  pair  are  used. 

The  data  required  for  a  subroutine  based  on  this  algorithm  are  q,  d,  S<> 
The  first  task  is  to  'elect  I  and  S*<s  Then  the  steps  given  in  the  following  box  are 
performed. 
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p”  Tp't;  j 


For  each  i  contained  in  S*0  perform  operations  1  through  5: 

1.  Replace  the  row  <Aj>in  A  with  <Aj>  -  <Aj>. 

2.  Similarly  replace  bi  with  bj  -  (bj  -d). 

3.  In  the  existing  vector,  b,  replace  {b}  with  {b}  -  —  { Aj} . 

cj 

4.  In  the  existing  matrix.  A,  replace  the  column  {Ai}  with  {Ai}  -  {Aj}. 

5.  Set  Ay  *  Cj,  Au  «  q,  An  *  q  and  bj  ■  d.  Place  zeros  in  the  remainin 
entries  of  column  I  and  row  I  of  A. 

For  each  i  and  j  contained  in  S*0  perform  the  following  operation: 

6.  Replace  the  component  Ajj  with  Ajj  +  . 


The  result  is  a  symmetric  system  whose  solution  automatically  satisfies  the 
multipoint  constraint. 


MULTIPLE  CONSTRAINTS 


If  several  constraints  exist  but  the  constraints  contain  variables  that  are  not 
common,  then  the  constraints  can  be  applied  sequentially  with  no  modifications  to 
the  proposed  algorithm. 

Frequently  it  is  necessary  to  apply  multiple  constraints  with  one  or  more 
variables  common  to  more  than  one  constraint.  With  the  proposed  approach  it  is 
possible  to  erase  a  previous  constraint  if  special  precautions  are  not  taken.  One 
method  to  circumvent  the  problem  is  to  scan  all  constraints  to  see  if,  for  each 
constraint,  a  variable  can  be  defined  which  does  not  appear  in  subsequent 
constraints.  Then  the  procedure  can  be  applied  sequentially  with  no 
modifications. 

If  it  proves  to  be  impossible  to  identify  a  unique  constraint  variable  for  a 
given  constraint,  then  an  alternative  approach  must  be  used.  When  this  situation 
arises,  the  constrain*  variables  must  be  eliminated  sequentially  using  previous 
constraint  equations  until  a  variable  appears  that  has  not  been  previously 
designated  a  constraint  variable.  Then  the  algorithm  can  be  applied  as  given. 


SIMPLE  EXAMPLE 

Consider  a  governing  matrix  and  a  force  vector  that  typically  result  from 
the  use  of  two-node  elements  in  one  dimension: 
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Suppose  the  following  two  constraints  are  applied: 

Xl  *  0  X6  -  X5  +  X4  -  X3  ■  1 

The  resulting  matrix,  force  vector  and  solution  vector  are: 
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Note  that  the  governing  matrix  remains  symmetric,  although  the  band  width  is 
larger,  and  the  solution  vector  satisfies  the  constraints. 


ONE-DIMENSIONAL  WAVE  PROPAGATION 
Formulation 

The  dynamic  equation  of  motion  for  a  bar  of  unit  cross-sectional  area  is 

pa+f*  *  Eu>xx  a  *  u,tt  (21) 


in  which  x  and  t  are  the  spatial  and  temporal  independent  variables,  respectively, 
the  dependent  variable  u  denotes  the  displacement,  E  is  Young’s  modulus,  which 

is  assumed  to  be  constant,  and  p  the  mass  density.  The  wave  speed  is  c  *  /E/p . 
The  force  per  unit  length  is  f*.  For  an  element  of  length  h,  the  element  stiffness 
and  mass  matrices  for  a  two-node  element  are: 
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(22) 


The  global  mass  and  stiffness  matrices  are  assembled  in  the  usual  manner  to 
obtain  the  spatially  discretized  version  of  (21)  as  follows: 

Ma  +  Ku  *  f  (23) 
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in  which  M  and  K  denote  the  system  mass  and  stiffness  matrices,  respectively,  and 
f  the  system  force  vector.  Suppose  (23)  is  discretized  in  time: 


MaB  +  Kn®  ■  f® 


in  which  n®  and  a®  denote  approximations  to  n  and  a,  respectively,  at  the  time  t® 
ns,  n«0, 1, 2,...  for  a  time  step  of  s.  Consider  a  general  trapezoidal  rule: 


*  v"  +  sfotja®*1  +  (l*ai)a“] 

■  u"  +  sja2vttfl  +  (l-a2)vl1] 


(25a) 

(25b) 


in  which  a  i  and  a  2  are  two  free  parameters  chosen  to  obtain  specific  integrators. 
When  (24)  and  (25)  are  combined,  the  resulting  system  integrator  for  updating 
the  velocity  vector  is 


in  which 


Av»+l 


A  =  M  +  ajct2S2K 

r®  *=  Bv®-sKu®  +aisf®+1  +  (1  -ai)sf" 
B  =  M-ai(l-a2)s2K 


(27a) 

(27b) 

(27c) 


The  updated  displacement  is  obtained  with  the  use  of  (25b).  The  algorithm  is 
implicit.  It  is  unconditionally  stable  if  the  integration  parameters  satisfy  the 
constraints  ai  =02^  0.5.  There  is  no  numerical  dissipation  if  01+02  =  1.  An 
explicit  integrator  is  obtained  by  summing  rows  of  [A]  to  obtain  a  diagonal  matrix. 

Suppose  constraints  on  u  are  converted  to  constraints  on  v.  Then  the 
multipoint  constraint  algorithm  need  only  be  applied  to  the  matrix  [A]  prior  to  an 
LU  decomposition,  and  the  implicit  integrator  is  applied  in  the  usual  manner. 

Static  Problem 

Consider  a  bar  of  length  L  =  1.28  fixed  at  the  left  end,  x  =  0.  At  the  right 
end,  suppose  a  force  of  magnitude  one  is  applied.  Also  suppose  that  the  third 
quarter  of  the  bar,  Of  £  x  &  0.75,  is  rigid.  If  the  bar  is  discretized  with  64 
uniform  elements,  then  the  rigid  zone  requires  16  constraint  equations,  namely, 
the  displacements  of  the  nodes  of  each  element  in  the  rigid  zone  are  equal.  For  E 
*  1,  the  solution  for  displacement  and  stress  in  the  bar  with  and  without  the 
inclusion  is  shown  in  Figs,  la  and  lb,  respectively.  As  expected,  the  strain  and 
consequently  the  stress  m  the  rigid  inclusion  is  predicted  to  be  zero  although  the 
stress  cannot  be  zero  if  equilibrium  is  to  be  satisfied. 


Dynamic  Problem 

Now  consider  the  bar  to  be  free  at  the  left  end  and  the  unit  force  at  the 
right  end  is  applied  instantaneously  at  t  «=  0.  The  solution  is  a  step  wave  which 
propagates  from  the  right  end  to  die  left  at  the  wave  velocity  of  c  =  1  for  p  *  1 
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and  E  «  1.  Solutions  are  shown  in  Fig.  2  for  a  1  *  <22  *  0.5  for  various 
configurations  and  s  «  0.5 h/c  where  h  denotes  the  length  of  the  element. 

Figure  2a  gives  the  numerical  solution  at  a  time  of  L/2c  for  the 
displacement  and  velocity  for  128  elements  with  no  rigid  inclusion.  Displayed  are 
typical  dispersive  features  such  as  overshoot  and  oscillations  which  are 
characteristic  of  the  numerical  algorithm  and  are  nonphysical.  A  simple  way  to 
handle  an  inclusion  is  to  use  an  artificially  high  stiffness.  Figure  2b  shows  the 
solution  for  the  same  mesh  and  wave  speed  but  with  a  Young’s  modulus  of  E*  « 
1024  used  to  simulate  a  rigid  inclusion  identical  to  that  used  for  the  static  problem. 
Figure  2c  shows  the  solution  to  the  same  rigid  inclusion  problem  but  with  a  time 
step  of  s  *  0.5h/c*  where  c*  is  the  wave  speed  based  on  E*.  This  small  time  step 
is  more  typical  of  those  required  with  an  explicit  integrator  and  represents  a 
severe  limitation  of  the  use  of  an  artificial  stiffness. 

Corresponding  solutions  obtained  with  the  use  of  the  constraint  algorithm 
are  given  in  Figs.  3a,  3b,  and  3c  for  64,  128  and  512  elements,  respectively.  When 
the  wave  transits  one-quarter  of  the  bar  and  impinges  on  the  rigid  segment, 
motion  is  immediately  realized  at  the  other  end  of  the  inclusion.  Therefore,  at  a 
time  of  L/2c  the  wave  has  actually  travelled  three-quarters  of  the  bar  rather  than 
one-half  the  length  which  is  the  distance  travelled  for  a  bar  with  no  inclusion.  For 
all  cases,  the  velocity  and  displacement  constraint  over  the  inclusion  is  shown. 
Some  of  the  oscillations  are  reduced  with  refinement  of  the  mesh  and  time  step. 

A  small  amount  of  dissipation  can  be  used  to  remove  spurious  oscillations. 
Results  are  given  in  Fig.  4  forai  =  a  2  =  0.52, 128  elements  and  s  =  0.5h/c.  Note 
that  the  wave  shapes  are  considerably  cleaner  for  both  the  homogeneous  bar  (Fig 
4a)  and  the  bar  with  a  rigid  inclusion  (Fig.  4b). 


EFFECT  ON  CONDITION  NUMBER 

For  a  constraint  equation,  if  the  largest  coefficient  in  absolute  value  is  set 
equal  to  one,  then  the  constraint  vector  has  been  scaled  so  that  the  infinity  norm 
is  one.  For  large  problems  it  may  be  desirable  to  scale  the  constraint  vectors  to 
minimize  the  condition  number  as  defined  by  the  ratio  of  maximum  to  minimum 
eigenvalues.  For  model  static  problems  such  as  that  whose  solution  is  shown  in 
Fig.  1,  all  constraints  were  scaled  by  a  factor  T) .  The  condition  number  was  then 
determined  for  various  values  of  “H  for  a  wide  range  in  the  number  of  elements 
and  in  the  size  of  a  rigid  inclusion.  Generally  speaking,  a  scale  factor  of  0.5/h  for  a 
uniform  mesh  provided  nearly  optimal  condition  numbers  for  a  large  number  of 
cases.  However,  it  is  not  obvious  what  the  scale  factors  should  be  for  a  set  of 
arbitraiy  constraints  without  actually  performing  numerical  investigations. 


CONCLUSIONS 

A  simple  algorithm  for  handling  one  or  more  multipoint  constraints  has 
been  presented.  Tne  dimension  of  the  governing  matrix  is  retained  as  is  symmetry 
and  positive  definiteness  when  these  properties  are  present.  The  change  in  band 
width  is  dictated  only  by  the  span  of  degrees  of  freedom  in  each  constraint 
equation.  The  application  to  both  static  and  dynamic  problems  is  straightforward 
and  should  prove  to  be  a  useful  computational  technique  for  the  numerous 
engineering  problems  in  which  constraints  are  important. 
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LIST  OF  FIGURES 


1.  Static  solution  for  displacement  and  stress  in  a  bar: 

(a)  without  a  rigid  inclusion, 

(b)  with  a  rigid  inclusion. 

2.  Plots  of  displacement  and  velocity  as  a  function  of  x  at  time  t  *  0-5L/c: 

(a)  no  inclusion  (s  »  OJh/c), 

(b)  the  simulation  of  the  inclusion  using  an  artificially  large  E  and  a  "large" 
time  step  of  s  *  OJh/c, 

(c)  the  simulation  of  the  inclusion  using  an  artificially  large  E  and  a  "small" 
time  step  of  s  *  0Jh/c\ 

3.  Plots  of  displacement  and  velocity  as  a  function  of  x  at  time  t  *  0.5L/c  obtained 

with  the  constraint  algorithm  (s  *  0.5h/c): 

(a)  32  elements, 

(b)  128  elements, 

(c)  256  elements. 

4.  Plots  of  displacement  and  velocity  as  a  function  of  x  at  time  t  *  0  JL/c  obtained 

with  a  small  amount  of  numerical  damping  (128  elements  and  s  *  0_5h/c): 

(a)  no  inclusion, 

(b)  constraint  algorithm  with  a  rigid  inclusion. 
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